Nothing
R/funcs.R
In hierNet: A Lasso for Hierarchical Interactions
Defines functions
hierNet.varimp
error.bars
plot.hierNet.cv
hierNet.cv
permute.rows
balanced.folds
twonorm
critf.logistic
predict.hierNet.logistic
ADMM4.Lagrangian
ggdescent.logistic
admm4.logistic
hierNet.logistic
ggdescent.c
Compute.phat.c
Compute.yhat.c
compute.full.interactions.c
compute.interactions.c
Objective.logistic
Objective
admm4
predict.hierNet.path
predict.hierNet
hierNet.path
print.hierNet.cv
print.hierNet.path
print.hierNet
hierNet
Documented in
compute.interactions.c
critf.logistic
hierNet
hierNet.cv
hierNet.logistic
hierNet.path
hierNet.varimp
Objective
Objective.logistic
plot.hierNet.cv
predict.hierNet
predict.hierNet.logistic
predict.hierNet.path
print.hierNet
print.hierNet.cv
print.hierNet.path
#' A Lasso for interactions
#'
#' One of the main functions in the hierNet package. Builds a regression model
#' with hierarchically constrained pairwise interactions. Required inputs are
#' an x matrix of features (the columns are the features) and a y vector of
#' values. Reasonably fast for moderate sized problems (100-200 variables). We
#' are currently working on an alternate algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param lam Regularization parameter (>0). L1 penalty param is \code{lam *
#' (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lam *
#' delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE.
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param center Should features be centered? Default TRUE; FALSE should rarely
#' be used. This option is available for special uses only
#' @param stand.main Should main effects be standardized? Default TRUE.
#' @param stand.int Should interactions be standardized? Default FALSE.
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter: threshold for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return \item{bp}{p-vector of estimated "positive part" main effect (p=#
#' features)} \item{bn}{p-vector of estimated "negative part" main effect;
#' overall main effect estimated coefficients are bp-bn } \item{th}{ Matrix of
#' estimated interaction coefficients, of dimension p by p. Note: when output
#' from hierNet is printed, th is symmetrized (set to (th+t(th))/2) for
#' simplicity.} \item{obj}{Value of objective function at minimum.}
#' \item{lam}{Value of lambda used} \item{type}{Type of model fit- "gaussian"
#' or "logistic" (binomial)} \item{mx}{ p-vector of column means of x}
#' \item{sx}{ p-vector of column standard deviations of x} \item{my}{ mean of
#' y} \item{mzz}{ column means of feature product matrix} \item{szz}{ column
#' standard deviations of feature product matrix} \item{call}{The call to
#' hierNet}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{predict.hierNet}, \link{hierNet.cv}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' # fit a single hierNet model
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet(x,y,lam=50)
#' print(fit)
#'
#' # try strong (rather than weak) hierarchy
#' fit=hierNet(x,y,lam=50, strong=TRUE)
#' print(fit)
#'
#' # a typical analysis including cross-validation
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#'
#' lamhat=fitcv$lamhat.1se
#' fit2=hierNet(x,y,lam=lamhat)
#' yhat=predict(fit2,x)
#'
#' @export hierNet
hierNet <- function(x, y, lam, delta=1e-8, strong=FALSE, diagonal=TRUE, aa=NULL, zz=NULL, center=TRUE, stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,
step=1, maxiter=2000, backtrack=0.2, tol=1e-5,
trace=0) {
# Main Hiernet function for fitting at a single parameter lambda.
# Note: L1 penalty terms have parameter lam.l1 = lambda * (1-delta)
# and L2 penalty has parameter lam.l2 = lambda * delta.
#
# stand.main and stand.int refer to scaling
stopifnot(nrow(x) == length(y), lam >= 0, delta >= 0, delta <= 1)
stopifnot(!is.null(step) && !is.null(maxiter))
if (strong) stopifnot(!is.null(niter))
stopifnot(class(y) == "numeric")
stopifnot(class(lam) == "numeric")
stopifnot(class(delta) == "numeric")
stopifnot(class(step) == "numeric", step > 0, maxiter > 0)
stopifnot(is.finite(x), is.finite(y), is.finite(lam), is.finite(delta))
this.call <- match.call()
if (!center) cat("WARNING: center=FALSE should almost never be used. This option is available for special uses only.", fill=TRUE)
# center and (maybe) scale variables
x <- scale(x, center=center, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale") # may be NULL
if (center) {
my <- mean(y)
y <- y - my
} else my <- NULL
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- compute.interactions.c(x, diagonal=diagonal)
}
if (is.matrix(zz)) {
zz <- scale(zz, center=center, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale") # may be NULL
zz <- as.numeric(zz)
} else {
mzz <- szz <- NULL
#cat("Provided zz is not a matrix, so it's assumed to be already centered.", fill=TRUE)
}
xnum <- as.numeric(x)
p <- ncol(x)
lam.l1 <- lam * (1 - delta)
lam.l2 <- lam * delta
if (strong) {
# strong hierarchy -- use ADMM4
if (is.null(rho)) rho <- as.numeric(nrow(x))
stopifnot(is.numeric(rho), is.finite(rho))
aa <- admm4(x, xnum, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, zz=zz,
rho=rho, niter=niter, aa=aa, sym.eps=sym.eps, # ADMM params
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol, # GG params
trace=trace)
# lack of symmetry in theta means that sometimes strong hierarchy will be (very slightly violated)
ii <- aa$bp + aa$bn == 0
# note aa$th[ii, ] = 0 since weak hierarchy holds for sure
if (sum(ii) > 0 & sum(ii) < p) {
thr <- max(abs(aa$th[!ii, ii]))
if (thr > 0) {
cat(" thr = ",thr, fill=TRUE)
if (thr > 1e-3)
warning("Had to change ADMM's 'th' by more than 0.001 to make strong hier hold! Increase niter (and/or rho). ")
aa$th[abs(aa$th) <= thr] <- 0
}
}
} else {
# weak hierarchy -- a single call to generalized gradient descent
if (is.null(aa)) {
aa <- list(th=matrix(0, p, p), bp=rep(0, p), bn=rep(0, p))
} else {
stopifnot(dim(aa$th) == c(p,p), length(aa$bp) == p, length(aa$bn) == p)
}
# this could be improved by not actually creating V...
V <- matrix(0, p, p)
rho <- 0
aa <- ggdescent.c(x=x, xnum=xnum, zz=zz, y=y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal,
rho=rho, V=V,
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol,
aa=aa, trace=trace)
}
aa$lam <- lam
aa$delta <- delta
aa$type <- "gaussian"
aa$diagonal <- diagonal
aa$strong <- strong
aa$obj <- Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=strong)
aa$step <- step
aa$maxiter <- maxiter
aa$backtrack <- backtrack
aa$tol <- tol
if (strong) {
# ADMM parameters:
aa$rho <- rho
aa$niter <- niter
aa$sym.eps <- sym.eps
}
aa$mx <- mx
aa$sx <- sx
aa$my <- my
aa$mzz <- mzz
aa$szz <- szz
aa$call <- this.call
class(aa) <- "hierNet"
return(aa)
}
print.hierNet <- function(x, ...) {
cat("Call:\n")
dput(x$call)
th=(x$th+t(x$th))/2
o2=colSums(th^2)!=0
b=x$bp-x$bn
o=b!=0
b=b[o]
if (any(o2)) {
# model has interactions
th=th[o,o2,drop=FALSE]
tight <- rowSums(abs(th)) >= x$bp[o] + x$bn[o] - 1e-9
tt <- rep("", length(tight))
tt[tight] <- "*"
mat=cbind(b,th)
mat=round(mat,4)
mat <- cbind(mat, tt)
cat("\n")
cat("Non-zero coefficients:",fill=T)
cat(" (Rows are predictors with nonzero main effects)",fill=T)
cat(" (1st column is main effect)", fill=T)
cat(" (Next columns are nonzero interactions of row predictor)", fill=T)
cat(" (Last column indicates whether hierarchy constraint is tight.)",fill=T)
cat("\n")
dimnames(mat)=list(as.character(which(o)),c("Main effect",as.character(which(o2)),"Tight?"))
print(mat, quote = FALSE)
} else {
mat <- matrix(round(b,4), length(b), 1)
cat("\n")
cat("Non-zero coefficients:",fill=T)
cat(" (No interactions in this model)",fill=T)
cat("\n")
dimnames(mat)=list(as.character(which(o)),"Main effect")
print(mat, quote = FALSE)
}
invisible()
}
print.hierNet.path <- function(x, ...) {
cat("Call:\n")
dput(x$call)
b=x$bp-x$bn
mat=cbind(round(x$lam,2),round(x$obj,2),colSums(b!=0),apply(x$th!=0,3,function(a) sum(diag(a)) + sum((a+t(a)!=0)[upper.tri(a)])))
dimnames(mat)=list(NULL,c("Lambda", "Objective", "Number of main effects","Number of interactions"))
cat("\n")
print(mat, quote = FALSE)
invisible()
}
print.hierNet.cv <- function(x, ...) {
cat("Call:\n")
dput(x$call)
mat=cbind(round(x$lamlist,2),x$nonzero,round(x$cv.err,2),round(x$cv.se,2))
dimnames(mat)=list(NULL,c("Lambda", "Number of nonzero","Mean CV error", "SE"))
cat("\n")
print(mat, quote = FALSE)
cat("\n")
cat(c("lamhat=",round(x$lamhat,2),"lamhat.1se=",round(x$lamhat.1se,2)),fill=T)
invisible()
}
#' Fit a path of hierNet models- lasso models with interactions
#'
#' One of the main functions in the hierNet package. Fits a path of hierNet
#' models over different values of the regularization parameter. Calls hierNet,
#' which builds a regression model with hierarchically constrained pairwise
#' interactions. Required inputs are an x matrix of features (the columns are
#' the features) and a y vector of values. Reasonably fast for moderate sized
#' problems (100-200 variables). We are currently working on an alternate
#' algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param lamlist Optional vector of values of lambda (the regularization
#' parameter). L1 penalty param is \code{lamdbda * (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lambda
#' * delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param minlam Optional minimum value for lambda
#' @param maxlam Optional maximum value for lambda
#' @param nlam Number of values of lambda to be tried
#' @param flmin Fraction of maxlam; minlam= flmin*maxlam. If computation is
#' slow, try increasing flmin to focus on the sparser part of the path
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param strong Flag specifying strong hierarchy (true) or weak hierarchy
#' (false). Default false
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing rho. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{bp}{p by nlam matrix of estimated "positive part" main effects
#' (p=#variables)}
#' \item{bn}{p by nlam matrix of estimated "negative part" main
#' effects}
#' \item{th}{p by p by nlam array of estimated interaction
#' coefficients}
#' \item{obj}{nlam values of objective function, one per lambda
#' value}
#' \item{lamlist}{Vector of values of lambda used}
#' \item{mx}{ p-vector
#' of column means of x}
#' \item{sx}{ p-vector of column standard deviations of
#' x}
#' \item{my}{ mean of y}
#' \item{mzz}{ column means of feature product matrix}
#' \item{szz}{ column standard deviations of feature product matrix}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{predict.hierNet}, \link{hierNet.cv}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' print(fit)
#'
#' @export hierNet.path
hierNet.path <- function(x, y, lamlist=NULL, delta=1e-8, minlam=NULL, maxlam=NULL, nlam=20, flmin=.01,
diagonal=TRUE, strong=FALSE, aa=NULL, zz=NULL,
stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG descent params
trace=0) {
# Main Hiernet function for fitting at a sequence of lambda values.
# Note: L1 penalty terms have parameter lam.l1 = lambda * (1-delta)
# and L2 penalty has parameter lam.l2 = lambda * delta.
#
# Always centers both x and zz (unless zz is provided in as.numeric form)
# stand.main and stand.int refer to whether main effects and interactions should have norm sqrt(n-1)
# center and (maybe) scale variables
this.call <- match.call()
x <- scale(x, center=TRUE, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale") # may be NULL
my <- mean(y)
y <- y - my
if (is.null(maxlam)) {
if (!is.null(minlam)) stop("Cannot have maxlam=NULL if minlam is non-null.")
# maxlam <- max(abs(t(x) %*% y)/colSums(x^2))
maxlam <- max(abs(t(x) %*% y))
# temp <- t(scale(t(x), center=FALSE, scale=1/y))
# temp2 <- apply(temp, 2, twonorm)
# maxlam <- max(max(temp2), maxlam)
minlam <- maxlam * flmin
}
if (is.null(minlam)) minlam <- maxlam * flmin
if (is.null(lamlist))
lamlist <- exp(seq(log(maxlam),log(minlam),length=nlam))
nlam <- length(lamlist)
if (is.null(zz))
zz <- compute.interactions.c(x, diagonal=diagonal)
else
stopifnot(is.matrix(zz))
# center and (maybe) scale zz
zz <- scale(zz, center=TRUE, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale") # may be NULL
zz <- as.numeric(zz)
p <- ncol(x)
cp2 <- choose(p, 2)
bp <- bn <- matrix(NA, nrow=p, ncol=nlam)
th <- array(NA, c(p, p, nlam))
obj <- rep(NA, nlam)
aa <- NULL
for (i in seq(nlam)) {
cat(c("i,lam=", i, round(lamlist[i],2)), fill=TRUE)
aa <- hierNet(x, y, lam=lamlist[i], delta=delta, strong=strong, diagonal=diagonal, aa=aa, zz=zz,
stand.main=FALSE, stand.int=FALSE, # have already standardized
rho=rho, niter=niter, sym.eps=sym.eps,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol, trace=trace)
bp[, i] <- aa$bp
bn[, i] <- aa$bn
th[, , i] <- aa$th
obj[i] <- aa$obj
}
dimnames(bp) <- dimnames(bn) <- list(as.character(1:p), NULL)
dimnames(th) <- list(as.character(1:p), as.character(1:p), NULL)
out <- list(bp=bp, bn=bn, th=th, obj=obj, lamlist=lamlist, delta=delta, mx=mx, sx=sx, mzz=mzz, szz=szz, my=my,
type="gaussian", diagonal=diagonal, strong=strong,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
call=this.call)
if (strong) {
# ADMM parameters:
out$rho <- rho
out$niter <- niter
out$sym.eps <- sym.eps
}
class(out) <- "hierNet.path"
out
}
#' Prediction function for hierNet and hierNet.logistic.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet" or "hiernet.logistic" function.
#'
#'
#' @param object The results of a call to the "hierNet" or "hierNet.path" or
#' function. The coefficients that are part of this object will be used for
#' making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one obseration per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Vector of predictions for each observation. For logistic
#' model, these are the estimated probabilities.}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet <- function(object, newx, newzz=NULL, ...) {
n <- nrow(newx)
if (is.null(object$sx))
newx <- scale(newx, center=object$mx, scale=FALSE)
else
newx <- scale(newx, center=object$mx, scale=object$sx)
if (is.null(newzz))
newzz <- compute.interactions.c(newx, diagonal=object$diagonal)
if (is.null(object$szz))
newzz <- scale(newzz, center=object$mzz, scale=FALSE)
else
newzz <- scale(newzz, center=object$mzz, scale=object$szz)
newzz <- as.numeric(newzz)
newx <- as.numeric(newx)
stopifnot(is.finite(newzz), is.finite(newx))
if (!("matrix" %in% class(object$bp)))
yhatt <- Compute.yhat.c(newx, newzz, object) + object$my
else {
nlam <- ncol(object$bp)
yhat <- matrix(NA, n, nlam)
# this could be made more efficient
for (i in seq(nlam)) {
bb <- list(bp=object$bp[, i], bn=object$bn[, i], th=object$th[, , i], diagonal=object$diagonal)
yhat[, i] <- Compute.yhat.c(newx, newzz, bb)
}
yhatt <- yhat + object$my
}
if (object$type == "logistic") {
# predict from hierNet.logistic object object
b0 <- object$b0
if(is.matrix(yhatt))
b0 <- matrix(b0, nrow=nrow(yhatt), ncol=ncol(yhatt), byrow=T)
yhatt <- b0 + yhatt
pr <- 1 / (1 + exp(-yhatt))
return(list(prob=pr, yhat=1*(pr>.5)))
}
return(yhatt)
}
#' Prediction function for hierNet.path and hierNet.logistic.path.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet.path" or "hiernet.logistic.path" functions.
#'
#'
#' @param object The results of a call to the "hierNet" or "hierNet.path" or
#' function. The coefficients that are part of this object will be used for
#' making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one obseration per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Matrix of predictions, one row per observation. For
#' logistic model, these are the estimated probabilities.}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet.path <- function(object, newx, newzz=NULL, ...){
predict.hierNet(object, newx, newzz, ...)
}
admm4 <- function(x, xnum, y, lam.l1, lam.l2, diagonal, zz=NULL, rho, niter, aa=NULL, sym.eps=1e-3, trace=1, ...) {
# Performs ADMM4.
# Note: xnum is the matrix x as a numeric. Both are passed to avoid having to call as.numeric too
# many times.
p <- ncol(x)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- as.numeric(compute.interactions.c(x, diagonal=diagonal))
}
else if ("matrix" %in% class(zz)) zz <- as.numeric(zz)
if (is.null(aa)) {
aa <- list(u=matrix(0, p, p),
th=matrix(0, p, p),
bp=rep(0, p),
bn=rep(0, p),
tt=matrix(0, p, p),
diagonal=diagonal)
} else {
stopifnot(diagonal == aa$diagonal)
}
if (is.null(aa$tt) || is.null(aa$u)) {
aa$tt <- 0.5 * (aa$th + t(aa$th))
aa$u <- matrix(0, p, p)
}
obj <- Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps)
ll <- NULL
for (i in seq(niter)) {
if (trace > 0) cat(i, " ")
ll <- c(ll, ADMM4.Lagrangian(aa, xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho))
V <- aa$u - rho * aa$tt
gg <- ggdescent.c(x, xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho, V, trace=trace-1, aa=aa, ...)
aa$th <- gg$th
aa$bp <- gg$bp
aa$bn <- gg$bn
aa$tt <- (aa$th + t(aa$th)) / 2 + (aa$u + t(aa$u)) / (2 * rho)
aa$u <- aa$u + rho * (aa$th - aa$tt)
obj <- c(obj, Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps))
if (trace > 0) cat(obj[i+1], fill=TRUE)
}
if (max(abs(aa$th-t(aa$th))) > sym.eps)
cat("Attention: th not symmetric within the desired sym.eps. Run ADMM for more iterations. And try increasing rho.")
aa$obj <- obj
aa$lagr <- ll
aa
}
Objective <- function(aa, x, y, lam.l1, lam.l2, xnum=NULL, zz=NULL, strong=TRUE, sym.eps=1e-3) {
# evaluates the NewYal objective at aa.
if (strong) {
if (max(aa$th-t(aa$th)) > sym.eps) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e-5)||any(aa$bn < -1e-5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Zero diagonal violated.", fill=TRUE)
return(Inf)
}
if (is.null(zz)) {
zz <- as.numeric(compute.interactions.c(x, diagonal=aa$diagonal))
}
if (is.null(xnum)) xnum <- as.numeric(x)
r <- y - Compute.yhat.c(xnum, zz, aa)
pen <- lam.l1 * sum(aa$bp + aa$bn) + lam.l1 * sum(abs(aa$th))/2 + lam.l1 * sum(abs(diag(aa$th)))/2
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
sum(r^2)/2 + pen
}
Objective.logistic <- function(aa, x, y, lam.l1, lam.l2, xnum=NULL, zz=NULL, strong=TRUE, sym.eps=1e-3) {
# evaluates the logistic hiernet objective at aa.
stopifnot(y %in% c(0,1))
stopifnot("diagonal" %in% names(aa))
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Diagonal of Theta is nonzero.", fill=TRUE)
return(Inf)
}
if (strong) {
if (max(aa$th-t(aa$th)) > sym.eps) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e5)||any(aa$bn < -1e5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (is.null(zz)) {
zz <- as.numeric(scale(compute.interactions.c(x, diagonal=aa$diagonal), center=TRUE, scale=FALSE))
}
if (is.matrix(zz)) zz <- as.numeric(zz)
if (is.null(xnum)) xnum <- as.numeric(x)
phat <- Compute.phat.c(xnum, zz, aa)
loss <- -sum(y*log(phat)) - sum((1-y)*log(1-phat))
pen <- lam.l1 * sum(aa$bp + aa$bn) + lam.l1 * sum(abs(aa$th))/2 + lam.l1 * sum(abs(diag(aa$th)))/2
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
loss + pen
}
compute.interactions.c <- function(x, diagonal=TRUE) {
# Returns (uncentered) n by cp2 matrix of interactions.
# The columns of zz are in standard order (11), 12,13,14,...,(22),23,...
# z's (jk)th column is x_j * x_k
n <- nrow(x)
p <- ncol(x)
cp2 <- p * (p - 1) / 2
if (diagonal) {
cp2 <- cp2 + p
out <- .C("ComputeInteractionsWithDiagWithIndices",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * cp2),
i1=as.integer(rep(0, cp2)),
i2=as.integer(rep(0, cp2)), PACKAGE="hierNet")
}
else {
out <- .C("ComputeInteractionsWithIndices",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * cp2),
i1=as.integer(rep(0, cp2)),
i2=as.integer(rep(0, cp2)), PACKAGE="hierNet")
}
z <- matrix(out$z, n, cp2)
rownames(z) <- rownames(x)
if (is.null(colnames(x))) {
colnames(z) <- paste(out$i1, out$i2, sep=":")
}
else {
colnames(z) <- paste(colnames(x)[out$i1], colnames(x)[out$i2], sep=":")
}
z
}
compute.full.interactions.c <- function(x) {
# Returns (uncentered) n by p^2 matrix of interactions.
# The columns of zz are in standard order 11,12,13,14,...,23,...
# z's (jk)th column is x_j * x_k
n <- nrow(x)
p <- ncol(x)
out <- .C("ComputeFullInteractions",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * p^2),
PACKAGE="hierNet")
matrix(out$z, n, p^2)
}
Compute.yhat.c <- function(xnum, zz, aa) {
# aa: list containing bp, bn, th, diagonal
# note: zz is the n by cp2 matrix, whereas z is the n by p^2 one.
p <- length(aa$bp)
n <- length(xnum) / p
stopifnot(n==round(n))
stopifnot("diagonal" %in% names(aa))
if (aa$diagonal) stopifnot(length(zz) == n * (choose(p,2) + p))
else stopifnot(length(zz) == n * choose(p,2))
out <- .C("compute_yhat_zz_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(aa$diagonal),
as.double(aa$th),
aa$bp,
aa$bn,
yhat=rep(0, n),
PACKAGE="hierNet")
out$yhat
}
Compute.phat.c <- function(xnum, zz, aa) {
# aa: list containing b0, bp, bn, th
# note: zz is the n by cp2 matrix, whereas z is the n by p^2 one.
stopifnot(c("b0","bp","bn","th","diagonal") %in% names(aa))
p <- length(aa$bp)
n <- length(xnum) / p
if (is.matrix(xnum)) xnum <- as.numeric(xnum)
stopifnot(n == round(n))
if (aa$diagonal) stopifnot(length(zz) == n * (choose(p,2) + p))
else stopifnot(length(zz) == n * choose(p,2))
#void compute_phat_zz_R(double *x, int *n, int *p, double *zz, int *diagonal,
# double *b0, double *th, double *bp, double *bn, double *phat) {
out <- .C("compute_phat_zz_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(aa$diagonal),
as.double(aa$b0),
as.double(aa$th),
aa$bp,
aa$bn,
phat=rep(0, n),
PACKAGE="hierNet")
out$phat
}
ggdescent.c <- function(x, xnum, zz, y, lam.l1, lam.l2, diagonal, rho, V, stepsize, backtrack=0.2, maxiter=100,
tol=1e-5, aa=NULL, trace=1) {
# See ADMM4 pdf for the problem this solves.
#
# x, xnum, zz, y: data (note: zz is a length n*cp2 vector, not a matrix) xnum is x as a vector
# lam.l1: l1-penalty parameter
# lam.l2: l2-penalty parameter
# rho: admm parameter
# V: see ADMM4 pdf
# stepsize: step size to start backtracking with
# backtrack: factor by which step is reduced on each backtrack.
# maxiter: number of generalized gradient steps to take.
# tol: stop gg descent if change in objective is below tol.
# aa: initial estimate of (th, bp, bn)
# trace: how verbose to be
#
# void ggdescent_R(double *x, int *n, int *p, double *zz, int *diagonal, double *y,
# double *lamL1, double*lamL2, double *rho, double *V, int *maxiter,
# double *curth, double *curbp, double *curbn,
# double *t, int *stepwindow, double *backtrack, double *tol, int *trace,
# double *th, double *bp, double *bn) {
n <- length(y)
p <- ncol(x)
stepwindow <- 10
if (is.null(aa)) aa <- list(th=matrix(0,p,p), bp=rep(0,p), bn=rep(0,p))
out <- .C("ggdescent_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(diagonal),
y,
as.double(lam.l1),
as.double(lam.l2),
as.double(rho),
as.double(V),
as.integer(maxiter),
as.double(aa$th),
aa$bp,
aa$bn,
stepsize,
as.integer(stepwindow),
backtrack,
tol,
as.integer(trace),
th=rep(0, p*p),
bp=rep(0, p),
bn=rep(0, p),
PACKAGE="hierNet")
list(bp=out$bp, bn=out$bn, th=matrix(out$th, p, p))
}
#' A logistic regression Lasso for interactions
#'
#' One of the main functions in the hierNet package. Builds a logistic
#' regression model with hierarchically constrained pairwise interactions.
#' Required inputs are an x matrix of features (the columns are the features)
#' and a y vector of values. Reasonably fast for moderate sized problems
#' (100-200 variables). We are currently working on a alternate algorithm for
#' large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, with values 0 or 1, where length(y)
#' equals nrow(x)
#' @param lam Regularization parameter (>0). L1 penalty param is \code{lam *
#' (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lam *
#' delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param center Should features be centered? Default TRUE; FALSE should rarely
#' be used. This option is available for special uses only
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{b0}{Intercept}
#' \item{bp}{p-vector of estimated "positive part"
#' main effect (p=#features)}
#' \item{bn}{p-vector of estimated "negative part"
#' main effect; overall main effect estimated coefficients are bp-bn }
#' \item{th}{Matrix of estimated interaction coefficients, of dimension p by p}
#' \item{obj}{Value of objective function at minimum.}
#' \item{lam}{Value of lambda used}
#' \item{type}{Type of model fit- "gaussian" or "logistic"
#' (binomial)}
#' \item{mx}{p-vector of column means of x}
#' \item{my}{Mean of y}
#' \item{sx}{p-vector of column standard deviations of x}
#' \item{mzz}{ column means of feature product matrix}
#' \item{call}{The call to hierNet}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{predict.hierNet.logistic},\link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic(x,y,lam=5)
#' print(fit)
#'
#' @export hierNet.logistic
hierNet.logistic <- function(x, y, lam, delta=1e-8, diagonal=TRUE, strong=FALSE, aa=NULL, zz=NULL, center=TRUE,
stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG descent params
trace=1) {
# Solves the logistic regression hiernet. Returns (b0, bp, bn, th)
this.call <- match.call()
n <- nrow(x)
p <- ncol(x)
stopifnot(y %in% c(0,1))
stopifnot(length(y) == n, lam >= 0, delta >= 0, delta <= 1)
stopifnot(!is.null(step) && !is.null(maxiter))
stopifnot(class(lam) == "numeric")
stopifnot(class(delta) == "numeric")
stopifnot(class(step) == "numeric", step > 0, maxiter > 0)
stopifnot(is.finite(x), is.finite(y), is.finite(lam), is.finite(delta))
lam.l1 <- lam * (1 - delta)
lam.l2 <- lam * delta
if (!center)
cat("WARNING: center=FALSE should almost never be used. This option is available for special uses only.", fill = TRUE)
x <- scale(x, center = center, scale = stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale")
if (is.null(aa)) aa <- list(b0=0, bp=rep(0, p), bn=rep(0, p), th=matrix(0, p, p), diagonal=diagonal)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- compute.interactions.c(x, diagonal=diagonal)
}
if (is.matrix(zz)) {
zz <- scale(zz, center=center, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale")
zz <- as.numeric(zz)
} else {
mzz <- szz <- NULL
#cat("Provided zz is not a matrix, so it's assumed to be already centered.", fill = TRUE)
}
xnum <- as.numeric(x)
if (strong) {
# strong hierarchy -- use ADMM4 (logistic regression version)
stopifnot(is.numeric(rho), is.finite(rho))
out <- admm4.logistic(x, xnum, y, lam.l1, lam.l2, diagonal=diagonal, zz=zz,
rho=rho, niter=niter, aa=aa, sym.eps=sym.eps, # ADMM params
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol, # GG params
trace=trace)
ii <- out$bp + out$bn == 0
# note out$th[ii, ] = 0 since weak hierarchy holds for sure
sumii <- sum(ii)
if (sumii > 0 && sumii < p) {
thr <- max(abs(out$th[!ii, ii]))
if (thr > 0) {
cat(" thr = ",thr, fill=TRUE)
if (thr > 1e-3)
warning("Had to change ADMM's 'th' by more than 0.001 to make strong hier hold! Increase niter (and/or rho). ")
aa$th[abs(aa$th) <= thr] <- 0
}
}
} else {
out <- ggdescent.logistic(xnum=xnum, zz=zz, y=y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho=0, V=matrix(0,p,p),
stepsize=step, backtrack=backtrack, maxiter=maxiter,
tol=tol, aa=aa, trace=trace)
}
out$call <- this.call
out$lam <- lam
out$delta <- delta
out$type <- "logistic"
out$diagonal <- diagonal
out$strong <- strong
if (strong) {
# ADMM parameters:
out$rho <- rho
out$niter <- niter
out$sym.eps <- sym.eps
}
out$step <- step
out$maxiter <- maxiter
out$backtrack <- backtrack
out$tol <- tol
out$obj <- critf.logistic(x, y, lam.l1, lam.l2, out$b0, out$bp, out$bn, out$th)
out$mx <- mx
out$my <- 0
out$sx <- sx
out$mzz <- mzz
class(out) <- "hierNet"
return(out)
}
admm4.logistic <- function(x, xnum, y, lam.l1, lam.l2, diagonal, zz=NULL, rho=10, niter, aa=NULL, sym.eps=1e-3, trace=1, ...) {
# Performs ADMM4 for logistic loss.
# Note: xnum is the matrix x as a numeric. Both are passed to avoid having to call as.numeric too
# many times.
p <- ncol(x)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- as.numeric(compute.interactions.c(x, diagonal=diagonal))
}
else if ("matrix" %in% class(zz)) zz <- as.numeric(zz)
if (is.null(aa)) {
aa <- list(u=matrix(0, p, p),
th=matrix(0, p, p),
bp=rep(0, p),
bn=rep(0, p),
tt=matrix(0, p, p),
diagonal=diagonal)
}
if (is.null(aa$tt) || is.null(aa$u)) {
aa$tt <- 0.5 * (aa$th + t(aa$th))
aa$u <- matrix(0, p, p)
}
obj <- Objective.logistic(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps)
for (i in seq(niter)) {
if (trace > 0) cat(i, " ")
V <- aa$u - rho * aa$tt
gg <- ggdescent.logistic(xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho, V, trace=trace-1, aa=aa, ...)
aa$th <- gg$th
aa$bp <- gg$bp
aa$bn <- gg$bn
aa$tt <- (aa$th + t(aa$th)) / 2 + (aa$u + t(aa$u)) / (2 * rho)
aa$u <- aa$u + rho * (aa$th - aa$tt)
obj <- c(obj, Objective.logistic(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps))
if (trace > 0) cat(obj[i+1], fill=TRUE)
}
if (max(abs(aa$th-t(aa$th))) > sym.eps)
cat("Attention: th not symmetric within the desired sym.eps. Run ADMM for more iterations. And try increasing rho.")
aa$obj <- obj
aa
}
ggdescent.logistic <- function(xnum, zz, y, lam.l1, lam.l2, diagonal, rho, V, stepsize, backtrack=0.2, maxiter=100,
tol=1e-5, aa=NULL, trace=1) {
# See ADMM4 pdf and logistic.pdf for the problem this solves.
#
# xnum, zz, y: data (note: zz is a length n*cp2 vector, not a matrix) xnum is x as a (n*p)-vector
# lam.l1: l1-penalty parameter
# lam.l2: l2-penalty parameter
# rho: admm parameter
# V: see ADMM4 pdf
# stepsize: step size to start backtracking with
# backtrack: factor by which step is reduced on each backtrack.
# maxiter: number of generalized gradient steps to take.
# tol: stop gg descent if change in objective is below tol.
# aa: initial estimate of (b0, th, bp, bn)
# trace: how verbose to be
#
#void ggdescent_logistic_R(double *x, int *n, int *p, double *zz, int * diagonal, double *y,
# double *lamL1, double *lamL2, double *rho, double *V, int *maxiter,
# double *curb0, double *curth, double *curbp, double *curbn,
# double *t, int *stepwindow, double *backtrack, double *tol, int *trace,
# double *b0, double *th, double *bp, double *bn) {
n <- length(y)
p <- length(xnum) / n
stopifnot(p == round(p))
if (diagonal) stopifnot(length(zz) == n * (choose(p,2)+p))
else stopifnot(length(zz) == n * choose(p,2))
stepwindow <- 10
if (is.null(aa)) aa <- list(b0=0, th=matrix(0,p,p), bp=rep(0,p), bn=rep(0,p))
out <- .C("ggdescent_logistic_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(diagonal),
as.double(y), # convert from integer to double
as.double(lam.l1),
as.double(lam.l2),
as.double(rho),
as.double(V),
as.integer(maxiter),
as.double(aa$b0),
as.double(aa$th),
aa$bp,
aa$bn,
as.double(stepsize),
as.integer(stepwindow),
as.double(backtrack),
as.double(tol),
as.integer(trace),
b0=as.double(0),
th=rep(0, p*p),
bp=rep(0, p),
bn=rep(0, p),
PACKAGE="hierNet")
list(b0=out$b0, bp=out$bp, bn=out$bn, th=matrix(out$th, p, p))
}
ADMM4.Lagrangian <- function(aa, xnum, zz, y, lam.l1, lam.l2, diagonal, rho) {
# aa: list with (th, bp, bn, tt, u)
# zz is a vector not a matrix
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Diagonal of Theta is nonzero.", fill=TRUE)
return(Inf)
}
if (max(aa$tt-t(aa$tt)) > 1e-8) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e-5)||any(aa$bn < -1e-5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-5)) {
cat("Zero diagonal violated.", fill=TRUE)
return(Inf)
}
V <- aa$u - rho * aa$tt
r <- y - Compute.yhat.c(xnum, zz, aa)
admm <- sum(aa$u*(aa$th-aa$tt)) + (rho/2) * sum((aa$th-aa$tt)^2)
#admm <- sum(V*aa$th) + (rho/2) * sum(aa$th^2) + (rho/2)*sum(aa$tt^2) - sum(aa$u*aa$tt)
pen <- lam.l1 * (sum(aa$bp + aa$bn) + sum(abs(aa$th))/2)
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
sum(r^2)/2 + pen + admm
}
#' Prediction function for hierNet.logistic.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet.logistic" function or "hierNet.logistic.path".
#'
#'
#' @param object The results of a call to the "hierNet.logistic" or
#' "hierNet.logistic.path" or function. The coefficients that are part of this
#' object will be used for making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one observation per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Matrix of predictions (probabilities), one row per
#' observation}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet.logistic}, \link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet.logistic(x,y,lam=5)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.logistic.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet.logistic <- function(object, newx, newzz=NULL, ...) {
predict.hierNet(object, newx, newzz, ...)
}
critf.logistic <- function(x, y, lam.l1, lam.l2, b0, bp, bn, th) {
yhat <- b0 + x %*% (bp - bn) + 0.5 * diag(x %*% th %*% t(x))
p <- 1 / (1 + exp(-yhat))
val <- -sum(y * log(p) + (1 - y) * log(1 - p))
val <- val + lam.l1 * sum(bp + bn) + lam.l1 * sum(abs(th))/2 + lam.l1 * sum(abs(diag(th)))/2
val <- val + lam.l2 * (sum(bp^2) + sum(bn^2) + sum(th^2))
return(val)
}
twonorm <- function(x) {sqrt(sum(x * x))}
#' Fit a path of logistic hierNet models- lasso models with interactions
#'
#' One of the main functions in the hierNet package. Fits a logistic path of
#' hierNet models over different values of the regularization parameter. Calls
#' hierNet.logistic, which builds a regression model with hierarchically
#' constrained pairwise interactions. Required inputs are an x matrix of
#' features (the columns are the features) and a y vector of values. Reasonably
#' fast for moderate sized problems (100-200 variables). We are currently
#' working on a alternate algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations equal to 0 or 1, where length(y) equals
#' nrow(x)
#' @param lamlist Optional vector of values of lambda (the regularization
#' parameter). L1 penalty param is \code{lambda * (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lambda
#' * delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param minlam Optional minimum value for lambda
#' @param maxlam Optional maximum value for lambda
#' @param flmin Fraction of maxlam; minlam= flmin*maxlam. If computation is
#' slow, try increasing flmin to focus on the sparser part of the path
#' @param nlam Number of values of lambda to be tried
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{bp}{p by nlam matrix of estimated "positive part" main effects
#' (p=#features)}
#' \item{bn}{p by nlam matrix of estimated "negative part" main effects}
#' \item{th}{p by p by nlam array of estimated interaction coefficients}
#' \item{obj}{nlam values of objective function, one per lambda value}
#' \item{lamlist}{Vector of values of lambda used}
#' \item{mx}{p-vector of column means of x}
#' \item{sx}{p-vector of column standard deviations of x}
#' \item{my}{mean of y}
#' \item{mzz}{column means of feature product matrix}
#' \item{szz}{column standard deviations of feature product matrix}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{predict.hierNet}, \link{hierNet.cv}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic.path(x,y)
#' print(fit)
#'
#' @export hierNet.logistic.path
hierNet.logistic.path <- function (x, y, lamlist=NULL, delta=1e-8, minlam=NULL, maxlam=NULL, flmin=.01, nlam=20,
diagonal=TRUE, strong=FALSE, aa=NULL,
zz=NULL, stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG params
trace=0) {
this.call=match.call()
stopifnot(y %in% c(0, 1))
x <- scale(x, center=TRUE, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale")
if (is.null(maxlam)) {
if (!is.null(minlam)) stop("Cannot have maxlam=NULL if minlam is non-null.")
maxlam <- max(abs(t(x) %*% y))
minlam <- maxlam * flmin
}
if (is.null(minlam)) minlam <- maxlam * flmin
if (is.null(lamlist))
lamlist <- exp(seq(log(maxlam), log(minlam), length=nlam))
nlam <- length(lamlist)
if (is.null(zz))
zz <- compute.interactions.c(x, diagonal=diagonal)
else stopifnot(is.matrix(zz))
zz <- scale(zz, center=TRUE, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale")
zz <- as.numeric(zz)
p <- ncol(x)
cp2 <- choose(p, 2)
b0 <- rep(NA, nlam)
bp <- bn <- matrix(NA, nrow=p, ncol=nlam)
th <- array(NA, c(p, p, nlam))
obj <- rep(NA, nlam)
aa <- NULL
for (i in seq(nlam)) {
cat(c("i,lam=", i, round(lamlist[i],2)), fill=TRUE)
aa <- hierNet.logistic(x, y, lam=lamlist[i], delta=delta, diagonal=diagonal, strong=strong,
aa=aa, zz=zz, stand.main=FALSE, stand.int=FALSE,
rho=rho, niter=niter, sym.eps=sym.eps,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
trace=trace)
b0[i] <- aa$b0
bp[, i] <- aa$bp
bn[, i] <- aa$bn
th[, , i] <- aa$th
obj[i] <- aa$obj
}
dimnames(bp) <- dimnames(bn) <- list(as.character(1:p), NULL)
dimnames(th) <- list(as.character(1:p), as.character(1:p), NULL)
out <- list(b0=b0, bp=bp, bn=bn, th=th, obj=obj, lamlist=lamlist, delta=delta,
mx=mx, my=0, sx=sx, mzz=mzz, szz=szz,
type="logistic", diagonal=diagonal, strong=strong,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
call=this.call)
if (strong) {
# ADMM parameters:
out$rho <- aa$rho
out$niter <- niter
out$sym.eps <- sym.eps
}
class(out) <- "hierNet.path"
out
}
balanced.folds <- function(y, nfolds=min(min(table(y)), 10)) {
totals <- table(y)
fmax <- max(totals)
nfolds <- min(nfolds, fmax)
# makes no sense to have more folds than the max class size
folds <- as.list(seq(nfolds))
yids <- split(seq(y), y)
# nice we to get the ids in a list, split by class
###Make a big matrix, with enough rows to get in all the folds per class
bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
for(i in seq(totals)) {
bigmat[seq(totals[i]), i] <- sample(yids[[i]])
}
smallmat <- matrix(bigmat, nrow = nfolds) # reshape the matrix
### Now do a clever sort to mix up the NAs
smallmat <- permute.rows(t(smallmat)) ### Now a clever unlisting
# the "clever" unlist doesn't work when there are no N
# apply(smallmat, 2, function(x)
# x[!is.na(x)])
res <-vector("list", nfolds)
for(j in 1:nfolds) {
jj <- !is.na(smallmat[, j])
res[[j]] <- smallmat[jj, j]
}
return(res)
}
permute.rows <-function(x) {
dd <- dim(x)
n <- dd[1]
p <- dd[2]
mm <- runif(length(x)) + rep(seq(n) * 10, rep(p, n))
matrix(t(x)[order(mm)], n, p, byrow = TRUE)
}
#' Cross-validation function for hierNet
#'
#' Uses cross-validation to estimate the regularization parameter for hierNet
#'
#'
#' @param fit Object returned from call to hierNet.path or
#' hierNet.logistic.path. All parameter settings will be taken from this
#' object.
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param nfolds Number of cross-validation folds
#' @param folds (Optional) user-supplied cross-validation folds. If provided,
#' nfolds is ignored.
#' @param trace Verbose output? 0=no, 1=yes
#' @return
#' \item{lamlist}{Vector of lambda values tried}
#' \item{cv.err}{Estimate of cross-validation error}
#' \item{cv.se}{Estimated standard error of cross-validation estimate}
#' \item{lamhat}{lambda value minimizing cv.err}
#' \item{lamhat.1se}{largest lambda value with cv.err less than or equal to
#' min(cv.err)+ SE }
#' \item{folds}{Indices of folds used in cross-validation}
#' \item{yhat}{n by nlam matrix of predicted values. Here, ith prediction is
#' based on training on all folds that do not include the ith data point.}
#' \item{nonzero}{Vector giving number of non-zero coefficients for each lambda
#' value}
#' \item{call}{The call to hierNet.cv}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{hierNet.path},
#' \link{hierNet.logistic},\link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#' plot(fitcv)
#'
#'
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#' plot(fitcv)
#'
#' @export hierNet.cv
hierNet.cv <- function(fit, x, y, nfolds=10, folds=NULL, trace=0) {
this.call <- match.call()
stopifnot(class(fit) == "hierNet.path")
if(fit$type=="gaussian"){errfun=function(y,yhat){(y-yhat)^2}}
if(fit$type=="logistic"){errfun=function(y,yhat){1*(y!=yhat)}}
n <- length(y)
if(is.null(folds)) {
folds <- split(sample(1:n), rep(1:nfolds, length = n))
}
else {
stopifnot(class(folds)=="list")
nfolds <- length(folds)
}
lamlist=fit$lamlist
# get whether fit was standardized based on fit$sx and fit$szz...
if (is.null(fit$mx)) stop("hierNet object was not centered. hierNet.cv has not been written for this (unusual) case.")
stand.main <- !is.null(fit$sx)
stand.int <- !is.null(fit$szz)
n.lamlist <- length(lamlist) ### Set up the data structures
size <- double(n.lamlist)
err2=matrix(NA,nrow=nfolds,ncol=length(lamlist))
for(ii in 1:nfolds) {
cat("Fold", ii, ":")
if(fit$type=="gaussian"){
a <- hierNet.path(x[-folds[[ii]],],y=y[-folds[[ii]]],
lamlist=lamlist, delta=fit$delta, diagonal=fit$diagonal, strong=fit$strong, trace=trace,
stand.main=stand.main, stand.int=stand.int,
rho=fit$rho, niter=fit$niter, sym.eps=fit$sym.eps, # ADMM parameters (which will be NULL if strong=F)
step=fit$step, maxiter=fit$maxiter, backtrack=fit$backtrack, tol=fit$tol) # GG descent params
yhatt=predict.hierNet(a,newx=x[folds[[ii]],])
}
if(fit$type=="logistic"){
a <- hierNet.logistic.path(x[-folds[[ii]],],y=y[-folds[[ii]]],
lamlist=lamlist, delta=fit$delta, diagonal=fit$diagonal, strong=fit$strong,
trace=trace, stand.main=stand.main, stand.int=stand.int,
rho=fit$rho, niter=fit$niter, sym.eps=fit$sym.eps, # ADMM parameters (which will be NULL if strong=F)
step=fit$step, maxiter=fit$maxiter, backtrack=fit$backtrack, tol=fit$tol) # GG descent params
yhatt=predict.hierNet.logistic(a,newx=x[folds[[ii]],])$yhat
}
temp=matrix(y[folds[[ii]]],nrow=length(folds[[ii]]),ncol=n.lamlist)
err2[ii,]=colMeans(errfun(yhatt,temp))
cat("\n")
}
errm=colMeans(err2)
errse=sqrt(apply(err2,2,var)/nfolds)
o=which.min(errm)
lamhat=lamlist[o]
oo=errm<= errm[o]+errse[o]
lamhat.1se=lamlist[oo & lamlist>=lamhat][1]
nonzero=colSums(fit$bp-fit$bn!=0) + apply(fit$th!=0, 3, function(a) sum(diag(a)) + sum((a+t(a)!=0)[upper.tri(a)]))
obj <- list(lamlist=lamlist, cv.err=errm,cv.se=errse,lamhat=lamhat, lamhat.1se=lamhat.1se,
nonzero=nonzero, folds=folds,
call = this.call)
class(obj) <- "hierNet.cv"
obj
}
plot.hierNet.cv <- function(x, ...) {
par(mar = c(5, 5, 5, 1))
yrang=range(c(x$cv.err-x$cv.se,x$cv.err+x$cv.se))
plot(log(x$lamlist), x$cv.err, xlab="log(lambda)",
ylab = "Cross-validation Error", type="n",ylim=yrang)
axis(3, at = log(x$lamlist), labels = paste(x$nonzero), srt = 90, adj = 0)
mtext("Number of features", 3, 4, cex = 1.2)
axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
error.bars(log(x$lamlist), x$cv.err - x$cv.se, x$cv.err + x$cv.se, width = 0.01, col = "darkgrey")
points(log(x$lamlist), x$cv.err, col=2, pch=19)
abline(v=log(x$lamhat), lty=3)
abline(v=log(x$lamhat.1se), lty=3)
invisible()
}
error.bars <-function(x, upper, lower, width = 0.02, ...) {
xlim <- range(x)
barw <- diff(xlim) * width
segments(x, upper, x, lower, ...)
segments(x - barw, upper, x + barw, upper, ...)
segments(x - barw, lower, x + barw, lower, ...)
range(upper, lower)
}
#' Variable importance for hierNet.
#'
#' (This is an experimental function.) Calculates a measure of the importance
#' of each variable.
#'
#'
#' @param fit The results of a call to the "hierNet"
#' @param x The training set feature matrix used in call produced "fit"
#' @param y The training set response vector used in call produced "fit"
#' @param ... additional arguments (not currently used)
#' @return Table of variable importance.
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
hierNet.varimp <- function(fit,x,y, ...) {
# NOTE: uses 0.5 cutoff for logistic case
lam=fit$lam
if(fit$type=="gaussian"){errfun=function(y,yhat){(y-yhat)^2}}
if(fit$type=="logistic"){
errfun=function(y,yhat){
term1=y*log(yhat);term1[yhat==0]=0
term2=(1-y)*log(1-yhat);term2[yhat==1]=0
val=-sum(term1+term2)
return(val)
}}
yhat=predict(fit,x)
rss=sum(errfun(y,yhat))
varsum=fit$bp-fit$bn+rowSums(abs(fit$th))
oo=which(abs(varsum)>1e-6)
imp=rss2=rep(NA,ncol(x))
for(j in oo){
cat(j)
fit0=fit;fit0$bp=fit$bp[-j];fit0$bn=fit$bn[-j];fit0$th=fit$th[-j,-j]
if(fit$type=="gaussian"){ fit2=hierNet(x[,-j],y,lam,delta=fit$delta,diagonal=fit$diagonal,aa=fit0)}
if(fit$type=="logistic"){ fit2=hierNet.logistic(x[,-j],y,lam,delta=fit$delta,diagonal=fit$diagonal,aa=fit0)}
yhat2=predict(fit2,x[,-j])
rss2[j]=sum(errfun(y,yhat2))
imp[j]=(rss2[j]-rss)/rss2[j]
}
imp[-oo]=0
res=cbind(1:ncol(x),round(imp,3))
ooo=order(-imp)
dimnames(res)=list(NULL,c("Predictor","Importance"))
cat("",fill=T)
return(res[ooo,])
}
Try the hierNet package in your browser
Any scripts or data that you put into this service are public.
hierNet documentation built on Sept. 10, 2026, 1:09 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions hierNet.varimp error.bars plot.hierNet.cv hierNet.cv permute.rows balanced.folds twonorm critf.logistic predict.hierNet.logistic ADMM4.Lagrangian ggdescent.logistic admm4.logistic hierNet.logistic ggdescent.c Compute.phat.c Compute.yhat.c compute.full.interactions.c compute.interactions.c Objective.logistic Objective admm4 predict.hierNet.path predict.hierNet hierNet.path print.hierNet.cv print.hierNet.path print.hierNet hierNet
Documented in compute.interactions.c critf.logistic hierNet hierNet.cv hierNet.logistic hierNet.path hierNet.varimp Objective Objective.logistic plot.hierNet.cv predict.hierNet predict.hierNet.logistic predict.hierNet.path print.hierNet print.hierNet.cv print.hierNet.path
#' A Lasso for interactions
#'
#' One of the main functions in the hierNet package. Builds a regression model
#' with hierarchically constrained pairwise interactions. Required inputs are
#' an x matrix of features (the columns are the features) and a y vector of
#' values. Reasonably fast for moderate sized problems (100-200 variables). We
#' are currently working on an alternate algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param lam Regularization parameter (>0). L1 penalty param is \code{lam *
#' (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lam *
#' delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE.
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param center Should features be centered? Default TRUE; FALSE should rarely
#' be used. This option is available for special uses only
#' @param stand.main Should main effects be standardized? Default TRUE.
#' @param stand.int Should interactions be standardized? Default FALSE.
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter: threshold for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return \item{bp}{p-vector of estimated "positive part" main effect (p=#
#' features)} \item{bn}{p-vector of estimated "negative part" main effect;
#' overall main effect estimated coefficients are bp-bn } \item{th}{ Matrix of
#' estimated interaction coefficients, of dimension p by p. Note: when output
#' from hierNet is printed, th is symmetrized (set to (th+t(th))/2) for
#' simplicity.} \item{obj}{Value of objective function at minimum.}
#' \item{lam}{Value of lambda used} \item{type}{Type of model fit- "gaussian"
#' or "logistic" (binomial)} \item{mx}{ p-vector of column means of x}
#' \item{sx}{ p-vector of column standard deviations of x} \item{my}{ mean of
#' y} \item{mzz}{ column means of feature product matrix} \item{szz}{ column
#' standard deviations of feature product matrix} \item{call}{The call to
#' hierNet}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{predict.hierNet}, \link{hierNet.cv}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' # fit a single hierNet model
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet(x,y,lam=50)
#' print(fit)
#'
#' # try strong (rather than weak) hierarchy
#' fit=hierNet(x,y,lam=50, strong=TRUE)
#' print(fit)
#'
#' # a typical analysis including cross-validation
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#'
#' lamhat=fitcv$lamhat.1se
#' fit2=hierNet(x,y,lam=lamhat)
#' yhat=predict(fit2,x)
#'
#' @export hierNet
hierNet <- function(x, y, lam, delta=1e-8, strong=FALSE, diagonal=TRUE, aa=NULL, zz=NULL, center=TRUE, stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,
step=1, maxiter=2000, backtrack=0.2, tol=1e-5,
trace=0) {
# Main Hiernet function for fitting at a single parameter lambda.
# Note: L1 penalty terms have parameter lam.l1 = lambda * (1-delta)
# and L2 penalty has parameter lam.l2 = lambda * delta.
#
# stand.main and stand.int refer to scaling
stopifnot(nrow(x) == length(y), lam >= 0, delta >= 0, delta <= 1)
stopifnot(!is.null(step) && !is.null(maxiter))
if (strong) stopifnot(!is.null(niter))
stopifnot(class(y) == "numeric")
stopifnot(class(lam) == "numeric")
stopifnot(class(delta) == "numeric")
stopifnot(class(step) == "numeric", step > 0, maxiter > 0)
stopifnot(is.finite(x), is.finite(y), is.finite(lam), is.finite(delta))
this.call <- match.call()
if (!center) cat("WARNING: center=FALSE should almost never be used. This option is available for special uses only.", fill=TRUE)
# center and (maybe) scale variables
x <- scale(x, center=center, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale") # may be NULL
if (center) {
my <- mean(y)
y <- y - my
} else my <- NULL
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- compute.interactions.c(x, diagonal=diagonal)
}
if (is.matrix(zz)) {
zz <- scale(zz, center=center, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale") # may be NULL
zz <- as.numeric(zz)
} else {
mzz <- szz <- NULL
#cat("Provided zz is not a matrix, so it's assumed to be already centered.", fill=TRUE)
}
xnum <- as.numeric(x)
p <- ncol(x)
lam.l1 <- lam * (1 - delta)
lam.l2 <- lam * delta
if (strong) {
# strong hierarchy -- use ADMM4
if (is.null(rho)) rho <- as.numeric(nrow(x))
stopifnot(is.numeric(rho), is.finite(rho))
aa <- admm4(x, xnum, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, zz=zz,
rho=rho, niter=niter, aa=aa, sym.eps=sym.eps, # ADMM params
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol, # GG params
trace=trace)
# lack of symmetry in theta means that sometimes strong hierarchy will be (very slightly violated)
ii <- aa$bp + aa$bn == 0
# note aa$th[ii, ] = 0 since weak hierarchy holds for sure
if (sum(ii) > 0 & sum(ii) < p) {
thr <- max(abs(aa$th[!ii, ii]))
if (thr > 0) {
cat(" thr = ",thr, fill=TRUE)
if (thr > 1e-3)
warning("Had to change ADMM's 'th' by more than 0.001 to make strong hier hold! Increase niter (and/or rho). ")
aa$th[abs(aa$th) <= thr] <- 0
}
}
} else {
# weak hierarchy -- a single call to generalized gradient descent
if (is.null(aa)) {
aa <- list(th=matrix(0, p, p), bp=rep(0, p), bn=rep(0, p))
} else {
stopifnot(dim(aa$th) == c(p,p), length(aa$bp) == p, length(aa$bn) == p)
}
# this could be improved by not actually creating V...
V <- matrix(0, p, p)
rho <- 0
aa <- ggdescent.c(x=x, xnum=xnum, zz=zz, y=y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal,
rho=rho, V=V,
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol,
aa=aa, trace=trace)
}
aa$lam <- lam
aa$delta <- delta
aa$type <- "gaussian"
aa$diagonal <- diagonal
aa$strong <- strong
aa$obj <- Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=strong)
aa$step <- step
aa$maxiter <- maxiter
aa$backtrack <- backtrack
aa$tol <- tol
if (strong) {
# ADMM parameters:
aa$rho <- rho
aa$niter <- niter
aa$sym.eps <- sym.eps
}
aa$mx <- mx
aa$sx <- sx
aa$my <- my
aa$mzz <- mzz
aa$szz <- szz
aa$call <- this.call
class(aa) <- "hierNet"
return(aa)
}
print.hierNet <- function(x, ...) {
cat("Call:\n")
dput(x$call)
th=(x$th+t(x$th))/2
o2=colSums(th^2)!=0
b=x$bp-x$bn
o=b!=0
b=b[o]
if (any(o2)) {
# model has interactions
th=th[o,o2,drop=FALSE]
tight <- rowSums(abs(th)) >= x$bp[o] + x$bn[o] - 1e-9
tt <- rep("", length(tight))
tt[tight] <- "*"
mat=cbind(b,th)
mat=round(mat,4)
mat <- cbind(mat, tt)
cat("\n")
cat("Non-zero coefficients:",fill=T)
cat(" (Rows are predictors with nonzero main effects)",fill=T)
cat(" (1st column is main effect)", fill=T)
cat(" (Next columns are nonzero interactions of row predictor)", fill=T)
cat(" (Last column indicates whether hierarchy constraint is tight.)",fill=T)
cat("\n")
dimnames(mat)=list(as.character(which(o)),c("Main effect",as.character(which(o2)),"Tight?"))
print(mat, quote = FALSE)
} else {
mat <- matrix(round(b,4), length(b), 1)
cat("\n")
cat("Non-zero coefficients:",fill=T)
cat(" (No interactions in this model)",fill=T)
cat("\n")
dimnames(mat)=list(as.character(which(o)),"Main effect")
print(mat, quote = FALSE)
}
invisible()
}
print.hierNet.path <- function(x, ...) {
cat("Call:\n")
dput(x$call)
b=x$bp-x$bn
mat=cbind(round(x$lam,2),round(x$obj,2),colSums(b!=0),apply(x$th!=0,3,function(a) sum(diag(a)) + sum((a+t(a)!=0)[upper.tri(a)])))
dimnames(mat)=list(NULL,c("Lambda", "Objective", "Number of main effects","Number of interactions"))
cat("\n")
print(mat, quote = FALSE)
invisible()
}
print.hierNet.cv <- function(x, ...) {
cat("Call:\n")
dput(x$call)
mat=cbind(round(x$lamlist,2),x$nonzero,round(x$cv.err,2),round(x$cv.se,2))
dimnames(mat)=list(NULL,c("Lambda", "Number of nonzero","Mean CV error", "SE"))
cat("\n")
print(mat, quote = FALSE)
cat("\n")
cat(c("lamhat=",round(x$lamhat,2),"lamhat.1se=",round(x$lamhat.1se,2)),fill=T)
invisible()
}
#' Fit a path of hierNet models- lasso models with interactions
#'
#' One of the main functions in the hierNet package. Fits a path of hierNet
#' models over different values of the regularization parameter. Calls hierNet,
#' which builds a regression model with hierarchically constrained pairwise
#' interactions. Required inputs are an x matrix of features (the columns are
#' the features) and a y vector of values. Reasonably fast for moderate sized
#' problems (100-200 variables). We are currently working on an alternate
#' algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param lamlist Optional vector of values of lambda (the regularization
#' parameter). L1 penalty param is \code{lamdbda * (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lambda
#' * delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param minlam Optional minimum value for lambda
#' @param maxlam Optional maximum value for lambda
#' @param nlam Number of values of lambda to be tried
#' @param flmin Fraction of maxlam; minlam= flmin*maxlam. If computation is
#' slow, try increasing flmin to focus on the sparser part of the path
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param strong Flag specifying strong hierarchy (true) or weak hierarchy
#' (false). Default false
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing rho. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{bp}{p by nlam matrix of estimated "positive part" main effects
#' (p=#variables)}
#' \item{bn}{p by nlam matrix of estimated "negative part" main
#' effects}
#' \item{th}{p by p by nlam array of estimated interaction
#' coefficients}
#' \item{obj}{nlam values of objective function, one per lambda
#' value}
#' \item{lamlist}{Vector of values of lambda used}
#' \item{mx}{ p-vector
#' of column means of x}
#' \item{sx}{ p-vector of column standard deviations of
#' x}
#' \item{my}{ mean of y}
#' \item{mzz}{ column means of feature product matrix}
#' \item{szz}{ column standard deviations of feature product matrix}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{predict.hierNet}, \link{hierNet.cv}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' print(fit)
#'
#' @export hierNet.path
hierNet.path <- function(x, y, lamlist=NULL, delta=1e-8, minlam=NULL, maxlam=NULL, nlam=20, flmin=.01,
diagonal=TRUE, strong=FALSE, aa=NULL, zz=NULL,
stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG descent params
trace=0) {
# Main Hiernet function for fitting at a sequence of lambda values.
# Note: L1 penalty terms have parameter lam.l1 = lambda * (1-delta)
# and L2 penalty has parameter lam.l2 = lambda * delta.
#
# Always centers both x and zz (unless zz is provided in as.numeric form)
# stand.main and stand.int refer to whether main effects and interactions should have norm sqrt(n-1)
# center and (maybe) scale variables
this.call <- match.call()
x <- scale(x, center=TRUE, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale") # may be NULL
my <- mean(y)
y <- y - my
if (is.null(maxlam)) {
if (!is.null(minlam)) stop("Cannot have maxlam=NULL if minlam is non-null.")
# maxlam <- max(abs(t(x) %*% y)/colSums(x^2))
maxlam <- max(abs(t(x) %*% y))
# temp <- t(scale(t(x), center=FALSE, scale=1/y))
# temp2 <- apply(temp, 2, twonorm)
# maxlam <- max(max(temp2), maxlam)
minlam <- maxlam * flmin
}
if (is.null(minlam)) minlam <- maxlam * flmin
if (is.null(lamlist))
lamlist <- exp(seq(log(maxlam),log(minlam),length=nlam))
nlam <- length(lamlist)
if (is.null(zz))
zz <- compute.interactions.c(x, diagonal=diagonal)
else
stopifnot(is.matrix(zz))
# center and (maybe) scale zz
zz <- scale(zz, center=TRUE, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale") # may be NULL
zz <- as.numeric(zz)
p <- ncol(x)
cp2 <- choose(p, 2)
bp <- bn <- matrix(NA, nrow=p, ncol=nlam)
th <- array(NA, c(p, p, nlam))
obj <- rep(NA, nlam)
aa <- NULL
for (i in seq(nlam)) {
cat(c("i,lam=", i, round(lamlist[i],2)), fill=TRUE)
aa <- hierNet(x, y, lam=lamlist[i], delta=delta, strong=strong, diagonal=diagonal, aa=aa, zz=zz,
stand.main=FALSE, stand.int=FALSE, # have already standardized
rho=rho, niter=niter, sym.eps=sym.eps,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol, trace=trace)
bp[, i] <- aa$bp
bn[, i] <- aa$bn
th[, , i] <- aa$th
obj[i] <- aa$obj
}
dimnames(bp) <- dimnames(bn) <- list(as.character(1:p), NULL)
dimnames(th) <- list(as.character(1:p), as.character(1:p), NULL)
out <- list(bp=bp, bn=bn, th=th, obj=obj, lamlist=lamlist, delta=delta, mx=mx, sx=sx, mzz=mzz, szz=szz, my=my,
type="gaussian", diagonal=diagonal, strong=strong,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
call=this.call)
if (strong) {
# ADMM parameters:
out$rho <- rho
out$niter <- niter
out$sym.eps <- sym.eps
}
class(out) <- "hierNet.path"
out
}
#' Prediction function for hierNet and hierNet.logistic.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet" or "hiernet.logistic" function.
#'
#'
#' @param object The results of a call to the "hierNet" or "hierNet.path" or
#' function. The coefficients that are part of this object will be used for
#' making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one obseration per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Vector of predictions for each observation. For logistic
#' model, these are the estimated probabilities.}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet <- function(object, newx, newzz=NULL, ...) {
n <- nrow(newx)
if (is.null(object$sx))
newx <- scale(newx, center=object$mx, scale=FALSE)
else
newx <- scale(newx, center=object$mx, scale=object$sx)
if (is.null(newzz))
newzz <- compute.interactions.c(newx, diagonal=object$diagonal)
if (is.null(object$szz))
newzz <- scale(newzz, center=object$mzz, scale=FALSE)
else
newzz <- scale(newzz, center=object$mzz, scale=object$szz)
newzz <- as.numeric(newzz)
newx <- as.numeric(newx)
stopifnot(is.finite(newzz), is.finite(newx))
if (!("matrix" %in% class(object$bp)))
yhatt <- Compute.yhat.c(newx, newzz, object) + object$my
else {
nlam <- ncol(object$bp)
yhat <- matrix(NA, n, nlam)
# this could be made more efficient
for (i in seq(nlam)) {
bb <- list(bp=object$bp[, i], bn=object$bn[, i], th=object$th[, , i], diagonal=object$diagonal)
yhat[, i] <- Compute.yhat.c(newx, newzz, bb)
}
yhatt <- yhat + object$my
}
if (object$type == "logistic") {
# predict from hierNet.logistic object object
b0 <- object$b0
if(is.matrix(yhatt))
b0 <- matrix(b0, nrow=nrow(yhatt), ncol=ncol(yhatt), byrow=T)
yhatt <- b0 + yhatt
pr <- 1 / (1 + exp(-yhatt))
return(list(prob=pr, yhat=1*(pr>.5)))
}
return(yhatt)
}
#' Prediction function for hierNet.path and hierNet.logistic.path.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet.path" or "hiernet.logistic.path" functions.
#'
#'
#' @param object The results of a call to the "hierNet" or "hierNet.path" or
#' function. The coefficients that are part of this object will be used for
#' making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one obseration per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Matrix of predictions, one row per observation. For
#' logistic model, these are the estimated probabilities.}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet.path <- function(object, newx, newzz=NULL, ...){
predict.hierNet(object, newx, newzz, ...)
}
admm4 <- function(x, xnum, y, lam.l1, lam.l2, diagonal, zz=NULL, rho, niter, aa=NULL, sym.eps=1e-3, trace=1, ...) {
# Performs ADMM4.
# Note: xnum is the matrix x as a numeric. Both are passed to avoid having to call as.numeric too
# many times.
p <- ncol(x)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- as.numeric(compute.interactions.c(x, diagonal=diagonal))
}
else if ("matrix" %in% class(zz)) zz <- as.numeric(zz)
if (is.null(aa)) {
aa <- list(u=matrix(0, p, p),
th=matrix(0, p, p),
bp=rep(0, p),
bn=rep(0, p),
tt=matrix(0, p, p),
diagonal=diagonal)
} else {
stopifnot(diagonal == aa$diagonal)
}
if (is.null(aa$tt) || is.null(aa$u)) {
aa$tt <- 0.5 * (aa$th + t(aa$th))
aa$u <- matrix(0, p, p)
}
obj <- Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps)
ll <- NULL
for (i in seq(niter)) {
if (trace > 0) cat(i, " ")
ll <- c(ll, ADMM4.Lagrangian(aa, xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho))
V <- aa$u - rho * aa$tt
gg <- ggdescent.c(x, xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho, V, trace=trace-1, aa=aa, ...)
aa$th <- gg$th
aa$bp <- gg$bp
aa$bn <- gg$bn
aa$tt <- (aa$th + t(aa$th)) / 2 + (aa$u + t(aa$u)) / (2 * rho)
aa$u <- aa$u + rho * (aa$th - aa$tt)
obj <- c(obj, Objective(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps))
if (trace > 0) cat(obj[i+1], fill=TRUE)
}
if (max(abs(aa$th-t(aa$th))) > sym.eps)
cat("Attention: th not symmetric within the desired sym.eps. Run ADMM for more iterations. And try increasing rho.")
aa$obj <- obj
aa$lagr <- ll
aa
}
Objective <- function(aa, x, y, lam.l1, lam.l2, xnum=NULL, zz=NULL, strong=TRUE, sym.eps=1e-3) {
# evaluates the NewYal objective at aa.
if (strong) {
if (max(aa$th-t(aa$th)) > sym.eps) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e-5)||any(aa$bn < -1e-5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Zero diagonal violated.", fill=TRUE)
return(Inf)
}
if (is.null(zz)) {
zz <- as.numeric(compute.interactions.c(x, diagonal=aa$diagonal))
}
if (is.null(xnum)) xnum <- as.numeric(x)
r <- y - Compute.yhat.c(xnum, zz, aa)
pen <- lam.l1 * sum(aa$bp + aa$bn) + lam.l1 * sum(abs(aa$th))/2 + lam.l1 * sum(abs(diag(aa$th)))/2
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
sum(r^2)/2 + pen
}
Objective.logistic <- function(aa, x, y, lam.l1, lam.l2, xnum=NULL, zz=NULL, strong=TRUE, sym.eps=1e-3) {
# evaluates the logistic hiernet objective at aa.
stopifnot(y %in% c(0,1))
stopifnot("diagonal" %in% names(aa))
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Diagonal of Theta is nonzero.", fill=TRUE)
return(Inf)
}
if (strong) {
if (max(aa$th-t(aa$th)) > sym.eps) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e5)||any(aa$bn < -1e5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (is.null(zz)) {
zz <- as.numeric(scale(compute.interactions.c(x, diagonal=aa$diagonal), center=TRUE, scale=FALSE))
}
if (is.matrix(zz)) zz <- as.numeric(zz)
if (is.null(xnum)) xnum <- as.numeric(x)
phat <- Compute.phat.c(xnum, zz, aa)
loss <- -sum(y*log(phat)) - sum((1-y)*log(1-phat))
pen <- lam.l1 * sum(aa$bp + aa$bn) + lam.l1 * sum(abs(aa$th))/2 + lam.l1 * sum(abs(diag(aa$th)))/2
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
loss + pen
}
compute.interactions.c <- function(x, diagonal=TRUE) {
# Returns (uncentered) n by cp2 matrix of interactions.
# The columns of zz are in standard order (11), 12,13,14,...,(22),23,...
# z's (jk)th column is x_j * x_k
n <- nrow(x)
p <- ncol(x)
cp2 <- p * (p - 1) / 2
if (diagonal) {
cp2 <- cp2 + p
out <- .C("ComputeInteractionsWithDiagWithIndices",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * cp2),
i1=as.integer(rep(0, cp2)),
i2=as.integer(rep(0, cp2)), PACKAGE="hierNet")
}
else {
out <- .C("ComputeInteractionsWithIndices",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * cp2),
i1=as.integer(rep(0, cp2)),
i2=as.integer(rep(0, cp2)), PACKAGE="hierNet")
}
z <- matrix(out$z, n, cp2)
rownames(z) <- rownames(x)
if (is.null(colnames(x))) {
colnames(z) <- paste(out$i1, out$i2, sep=":")
}
else {
colnames(z) <- paste(colnames(x)[out$i1], colnames(x)[out$i2], sep=":")
}
z
}
compute.full.interactions.c <- function(x) {
# Returns (uncentered) n by p^2 matrix of interactions.
# The columns of zz are in standard order 11,12,13,14,...,23,...
# z's (jk)th column is x_j * x_k
n <- nrow(x)
p <- ncol(x)
out <- .C("ComputeFullInteractions",
as.double(x),
as.integer(n),
as.integer(p),
z=rep(0, n * p^2),
PACKAGE="hierNet")
matrix(out$z, n, p^2)
}
Compute.yhat.c <- function(xnum, zz, aa) {
# aa: list containing bp, bn, th, diagonal
# note: zz is the n by cp2 matrix, whereas z is the n by p^2 one.
p <- length(aa$bp)
n <- length(xnum) / p
stopifnot(n==round(n))
stopifnot("diagonal" %in% names(aa))
if (aa$diagonal) stopifnot(length(zz) == n * (choose(p,2) + p))
else stopifnot(length(zz) == n * choose(p,2))
out <- .C("compute_yhat_zz_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(aa$diagonal),
as.double(aa$th),
aa$bp,
aa$bn,
yhat=rep(0, n),
PACKAGE="hierNet")
out$yhat
}
Compute.phat.c <- function(xnum, zz, aa) {
# aa: list containing b0, bp, bn, th
# note: zz is the n by cp2 matrix, whereas z is the n by p^2 one.
stopifnot(c("b0","bp","bn","th","diagonal") %in% names(aa))
p <- length(aa$bp)
n <- length(xnum) / p
if (is.matrix(xnum)) xnum <- as.numeric(xnum)
stopifnot(n == round(n))
if (aa$diagonal) stopifnot(length(zz) == n * (choose(p,2) + p))
else stopifnot(length(zz) == n * choose(p,2))
#void compute_phat_zz_R(double *x, int *n, int *p, double *zz, int *diagonal,
# double *b0, double *th, double *bp, double *bn, double *phat) {
out <- .C("compute_phat_zz_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(aa$diagonal),
as.double(aa$b0),
as.double(aa$th),
aa$bp,
aa$bn,
phat=rep(0, n),
PACKAGE="hierNet")
out$phat
}
ggdescent.c <- function(x, xnum, zz, y, lam.l1, lam.l2, diagonal, rho, V, stepsize, backtrack=0.2, maxiter=100,
tol=1e-5, aa=NULL, trace=1) {
# See ADMM4 pdf for the problem this solves.
#
# x, xnum, zz, y: data (note: zz is a length n*cp2 vector, not a matrix) xnum is x as a vector
# lam.l1: l1-penalty parameter
# lam.l2: l2-penalty parameter
# rho: admm parameter
# V: see ADMM4 pdf
# stepsize: step size to start backtracking with
# backtrack: factor by which step is reduced on each backtrack.
# maxiter: number of generalized gradient steps to take.
# tol: stop gg descent if change in objective is below tol.
# aa: initial estimate of (th, bp, bn)
# trace: how verbose to be
#
# void ggdescent_R(double *x, int *n, int *p, double *zz, int *diagonal, double *y,
# double *lamL1, double*lamL2, double *rho, double *V, int *maxiter,
# double *curth, double *curbp, double *curbn,
# double *t, int *stepwindow, double *backtrack, double *tol, int *trace,
# double *th, double *bp, double *bn) {
n <- length(y)
p <- ncol(x)
stepwindow <- 10
if (is.null(aa)) aa <- list(th=matrix(0,p,p), bp=rep(0,p), bn=rep(0,p))
out <- .C("ggdescent_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(diagonal),
y,
as.double(lam.l1),
as.double(lam.l2),
as.double(rho),
as.double(V),
as.integer(maxiter),
as.double(aa$th),
aa$bp,
aa$bn,
stepsize,
as.integer(stepwindow),
backtrack,
tol,
as.integer(trace),
th=rep(0, p*p),
bp=rep(0, p),
bn=rep(0, p),
PACKAGE="hierNet")
list(bp=out$bp, bn=out$bn, th=matrix(out$th, p, p))
}
#' A logistic regression Lasso for interactions
#'
#' One of the main functions in the hierNet package. Builds a logistic
#' regression model with hierarchically constrained pairwise interactions.
#' Required inputs are an x matrix of features (the columns are the features)
#' and a y vector of values. Reasonably fast for moderate sized problems
#' (100-200 variables). We are currently working on a alternate algorithm for
#' large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, with values 0 or 1, where length(y)
#' equals nrow(x)
#' @param lam Regularization parameter (>0). L1 penalty param is \code{lam *
#' (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lam *
#' delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param center Should features be centered? Default TRUE; FALSE should rarely
#' be used. This option is available for special uses only
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{b0}{Intercept}
#' \item{bp}{p-vector of estimated "positive part"
#' main effect (p=#features)}
#' \item{bn}{p-vector of estimated "negative part"
#' main effect; overall main effect estimated coefficients are bp-bn }
#' \item{th}{Matrix of estimated interaction coefficients, of dimension p by p}
#' \item{obj}{Value of objective function at minimum.}
#' \item{lam}{Value of lambda used}
#' \item{type}{Type of model fit- "gaussian" or "logistic"
#' (binomial)}
#' \item{mx}{p-vector of column means of x}
#' \item{my}{Mean of y}
#' \item{sx}{p-vector of column standard deviations of x}
#' \item{mzz}{ column means of feature product matrix}
#' \item{call}{The call to hierNet}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{predict.hierNet.logistic},\link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic(x,y,lam=5)
#' print(fit)
#'
#' @export hierNet.logistic
hierNet.logistic <- function(x, y, lam, delta=1e-8, diagonal=TRUE, strong=FALSE, aa=NULL, zz=NULL, center=TRUE,
stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG descent params
trace=1) {
# Solves the logistic regression hiernet. Returns (b0, bp, bn, th)
this.call <- match.call()
n <- nrow(x)
p <- ncol(x)
stopifnot(y %in% c(0,1))
stopifnot(length(y) == n, lam >= 0, delta >= 0, delta <= 1)
stopifnot(!is.null(step) && !is.null(maxiter))
stopifnot(class(lam) == "numeric")
stopifnot(class(delta) == "numeric")
stopifnot(class(step) == "numeric", step > 0, maxiter > 0)
stopifnot(is.finite(x), is.finite(y), is.finite(lam), is.finite(delta))
lam.l1 <- lam * (1 - delta)
lam.l2 <- lam * delta
if (!center)
cat("WARNING: center=FALSE should almost never be used. This option is available for special uses only.", fill = TRUE)
x <- scale(x, center = center, scale = stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale")
if (is.null(aa)) aa <- list(b0=0, bp=rep(0, p), bn=rep(0, p), th=matrix(0, p, p), diagonal=diagonal)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- compute.interactions.c(x, diagonal=diagonal)
}
if (is.matrix(zz)) {
zz <- scale(zz, center=center, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale")
zz <- as.numeric(zz)
} else {
mzz <- szz <- NULL
#cat("Provided zz is not a matrix, so it's assumed to be already centered.", fill = TRUE)
}
xnum <- as.numeric(x)
if (strong) {
# strong hierarchy -- use ADMM4 (logistic regression version)
stopifnot(is.numeric(rho), is.finite(rho))
out <- admm4.logistic(x, xnum, y, lam.l1, lam.l2, diagonal=diagonal, zz=zz,
rho=rho, niter=niter, aa=aa, sym.eps=sym.eps, # ADMM params
stepsize=step, backtrack=backtrack, maxiter=maxiter, tol=tol, # GG params
trace=trace)
ii <- out$bp + out$bn == 0
# note out$th[ii, ] = 0 since weak hierarchy holds for sure
sumii <- sum(ii)
if (sumii > 0 && sumii < p) {
thr <- max(abs(out$th[!ii, ii]))
if (thr > 0) {
cat(" thr = ",thr, fill=TRUE)
if (thr > 1e-3)
warning("Had to change ADMM's 'th' by more than 0.001 to make strong hier hold! Increase niter (and/or rho). ")
aa$th[abs(aa$th) <= thr] <- 0
}
}
} else {
out <- ggdescent.logistic(xnum=xnum, zz=zz, y=y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho=0, V=matrix(0,p,p),
stepsize=step, backtrack=backtrack, maxiter=maxiter,
tol=tol, aa=aa, trace=trace)
}
out$call <- this.call
out$lam <- lam
out$delta <- delta
out$type <- "logistic"
out$diagonal <- diagonal
out$strong <- strong
if (strong) {
# ADMM parameters:
out$rho <- rho
out$niter <- niter
out$sym.eps <- sym.eps
}
out$step <- step
out$maxiter <- maxiter
out$backtrack <- backtrack
out$tol <- tol
out$obj <- critf.logistic(x, y, lam.l1, lam.l2, out$b0, out$bp, out$bn, out$th)
out$mx <- mx
out$my <- 0
out$sx <- sx
out$mzz <- mzz
class(out) <- "hierNet"
return(out)
}
admm4.logistic <- function(x, xnum, y, lam.l1, lam.l2, diagonal, zz=NULL, rho=10, niter, aa=NULL, sym.eps=1e-3, trace=1, ...) {
# Performs ADMM4 for logistic loss.
# Note: xnum is the matrix x as a numeric. Both are passed to avoid having to call as.numeric too
# many times.
p <- ncol(x)
if (is.null(zz)) {
if (trace > 0) cat("Computing zz...", fill=TRUE)
zz <- as.numeric(compute.interactions.c(x, diagonal=diagonal))
}
else if ("matrix" %in% class(zz)) zz <- as.numeric(zz)
if (is.null(aa)) {
aa <- list(u=matrix(0, p, p),
th=matrix(0, p, p),
bp=rep(0, p),
bn=rep(0, p),
tt=matrix(0, p, p),
diagonal=diagonal)
}
if (is.null(aa$tt) || is.null(aa$u)) {
aa$tt <- 0.5 * (aa$th + t(aa$th))
aa$u <- matrix(0, p, p)
}
obj <- Objective.logistic(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps)
for (i in seq(niter)) {
if (trace > 0) cat(i, " ")
V <- aa$u - rho * aa$tt
gg <- ggdescent.logistic(xnum, zz, y, lam.l1=lam.l1, lam.l2=lam.l2, diagonal=diagonal, rho, V, trace=trace-1, aa=aa, ...)
aa$th <- gg$th
aa$bp <- gg$bp
aa$bn <- gg$bn
aa$tt <- (aa$th + t(aa$th)) / 2 + (aa$u + t(aa$u)) / (2 * rho)
aa$u <- aa$u + rho * (aa$th - aa$tt)
obj <- c(obj, Objective.logistic(aa=aa, x=x, y=y, lam.l1=lam.l1, lam.l2=lam.l2, xnum=xnum, zz=zz, strong=TRUE, sym.eps=sym.eps))
if (trace > 0) cat(obj[i+1], fill=TRUE)
}
if (max(abs(aa$th-t(aa$th))) > sym.eps)
cat("Attention: th not symmetric within the desired sym.eps. Run ADMM for more iterations. And try increasing rho.")
aa$obj <- obj
aa
}
ggdescent.logistic <- function(xnum, zz, y, lam.l1, lam.l2, diagonal, rho, V, stepsize, backtrack=0.2, maxiter=100,
tol=1e-5, aa=NULL, trace=1) {
# See ADMM4 pdf and logistic.pdf for the problem this solves.
#
# xnum, zz, y: data (note: zz is a length n*cp2 vector, not a matrix) xnum is x as a (n*p)-vector
# lam.l1: l1-penalty parameter
# lam.l2: l2-penalty parameter
# rho: admm parameter
# V: see ADMM4 pdf
# stepsize: step size to start backtracking with
# backtrack: factor by which step is reduced on each backtrack.
# maxiter: number of generalized gradient steps to take.
# tol: stop gg descent if change in objective is below tol.
# aa: initial estimate of (b0, th, bp, bn)
# trace: how verbose to be
#
#void ggdescent_logistic_R(double *x, int *n, int *p, double *zz, int * diagonal, double *y,
# double *lamL1, double *lamL2, double *rho, double *V, int *maxiter,
# double *curb0, double *curth, double *curbp, double *curbn,
# double *t, int *stepwindow, double *backtrack, double *tol, int *trace,
# double *b0, double *th, double *bp, double *bn) {
n <- length(y)
p <- length(xnum) / n
stopifnot(p == round(p))
if (diagonal) stopifnot(length(zz) == n * (choose(p,2)+p))
else stopifnot(length(zz) == n * choose(p,2))
stepwindow <- 10
if (is.null(aa)) aa <- list(b0=0, th=matrix(0,p,p), bp=rep(0,p), bn=rep(0,p))
out <- .C("ggdescent_logistic_R",
xnum,
as.integer(n),
as.integer(p),
zz,
as.integer(diagonal),
as.double(y), # convert from integer to double
as.double(lam.l1),
as.double(lam.l2),
as.double(rho),
as.double(V),
as.integer(maxiter),
as.double(aa$b0),
as.double(aa$th),
aa$bp,
aa$bn,
as.double(stepsize),
as.integer(stepwindow),
as.double(backtrack),
as.double(tol),
as.integer(trace),
b0=as.double(0),
th=rep(0, p*p),
bp=rep(0, p),
bn=rep(0, p),
PACKAGE="hierNet")
list(b0=out$b0, bp=out$bp, bn=out$bn, th=matrix(out$th, p, p))
}
ADMM4.Lagrangian <- function(aa, xnum, zz, y, lam.l1, lam.l2, diagonal, rho) {
# aa: list with (th, bp, bn, tt, u)
# zz is a vector not a matrix
if (aa$diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-8)) {
cat("Diagonal of Theta is nonzero.", fill=TRUE)
return(Inf)
}
if (max(aa$tt-t(aa$tt)) > 1e-8) {
cat("Theta is not symmetric.", fill=TRUE)
return(Inf)
}
if (any(rowSums(abs(aa$th)) > aa$bp + aa$bn + 1e-5)) {
cat("hierarchy violated.", fill=TRUE)
return(Inf)
}
if (any(aa$bp < -1e-5)||any(aa$bn < -1e-5)) {
cat("Non-negative of bp or bn violated.", fill=TRUE)
return(Inf)
}
if (diagonal == FALSE)
if (any(abs(diag(aa$th)) > 1e-5)) {
cat("Zero diagonal violated.", fill=TRUE)
return(Inf)
}
V <- aa$u - rho * aa$tt
r <- y - Compute.yhat.c(xnum, zz, aa)
admm <- sum(aa$u*(aa$th-aa$tt)) + (rho/2) * sum((aa$th-aa$tt)^2)
#admm <- sum(V*aa$th) + (rho/2) * sum(aa$th^2) + (rho/2)*sum(aa$tt^2) - sum(aa$u*aa$tt)
pen <- lam.l1 * (sum(aa$bp + aa$bn) + sum(abs(aa$th))/2)
pen <- pen + lam.l2 * (sum(aa$bp^2) + sum(aa$bn^2) + sum(aa$th^2))
sum(r^2)/2 + pen + admm
}
#' Prediction function for hierNet.logistic.
#'
#' A function to perform prediction, using an x matrix and the output of the
#' "hierNet.logistic" function or "hierNet.logistic.path".
#'
#'
#' @param object The results of a call to the "hierNet.logistic" or
#' "hierNet.logistic.path" or function. The coefficients that are part of this
#' object will be used for making predictions.
#' @param newx The new x at which predictions should be made. Can be a vector
#' or a matrix (one observation per row).
#' @param newzz Optional matrix of products of columns of newx, computed by
#' compute.interactions.c
#' @param ... additional arguments (not currently used)
#' @return \item{yhat}{Matrix of predictions (probabilities), one row per
#' observation}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet.logistic}, \link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet.logistic(x,y,lam=5)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.logistic.path(x,y)
#' yhat=predict(fit,newx)
#'
predict.hierNet.logistic <- function(object, newx, newzz=NULL, ...) {
predict.hierNet(object, newx, newzz, ...)
}
critf.logistic <- function(x, y, lam.l1, lam.l2, b0, bp, bn, th) {
yhat <- b0 + x %*% (bp - bn) + 0.5 * diag(x %*% th %*% t(x))
p <- 1 / (1 + exp(-yhat))
val <- -sum(y * log(p) + (1 - y) * log(1 - p))
val <- val + lam.l1 * sum(bp + bn) + lam.l1 * sum(abs(th))/2 + lam.l1 * sum(abs(diag(th)))/2
val <- val + lam.l2 * (sum(bp^2) + sum(bn^2) + sum(th^2))
return(val)
}
twonorm <- function(x) {sqrt(sum(x * x))}
#' Fit a path of logistic hierNet models- lasso models with interactions
#'
#' One of the main functions in the hierNet package. Fits a logistic path of
#' hierNet models over different values of the regularization parameter. Calls
#' hierNet.logistic, which builds a regression model with hierarchically
#' constrained pairwise interactions. Required inputs are an x matrix of
#' features (the columns are the features) and a y vector of values. Reasonably
#' fast for moderate sized problems (100-200 variables). We are currently
#' working on a alternate algorithm for large scale problems.
#'
#'
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations equal to 0 or 1, where length(y) equals
#' nrow(x)
#' @param lamlist Optional vector of values of lambda (the regularization
#' parameter). L1 penalty param is \code{lambda * (1-delta)}.
#' @param delta Elastic Net parameter. Squared L2 penalty param is \code{lambda
#' * delta}. Not a tuning parameter: Think of as fixed and small. Default 1e-8.
#' @param minlam Optional minimum value for lambda
#' @param maxlam Optional maximum value for lambda
#' @param flmin Fraction of maxlam; minlam= flmin*maxlam. If computation is
#' slow, try increasing flmin to focus on the sparser part of the path
#' @param nlam Number of values of lambda to be tried
#' @param diagonal Flag specifying whether to include "pure" quadratic terms,
#' th_jjX_j^2, in the model. Default TRUE.
#' @param stand.main Should main effects be standardized? Default TRUE
#' @param stand.int Should interactions be standardized? Default FALSE
#' @param strong Flag specifying strong hierarchy (TRUE) or weak hierarchy
#' (FALSE). Default FALSE
#' @param aa An *optional* argument, a list with results from a previous call
#' @param zz An *optional* argument, a matrix whose columns are products of
#' features, computed by the function compute.interactions.c
#' @param rho ADMM parameter: tuning parameter (>0) for ADMM. If there are
#' convergence problems, try decreasing \code{rho}. Default n.
#' @param niter ADMM parameter: number of iterations
#' @param sym.eps ADMM parameter Thresholding for symmetrizing with strong=TRUE
#' @param step Stepsize for generalized gradient descent
#' @param maxiter Maximum number of iterations for generalized gradient descent
#' @param backtrack Backtrack parameter for generalized gradient descent
#' @param tol Error tolerance parameter for generalized gradient descent
#' @param trace Output option; trace=1 gives verbose output
#' @return
#' \item{bp}{p by nlam matrix of estimated "positive part" main effects
#' (p=#features)}
#' \item{bn}{p by nlam matrix of estimated "negative part" main effects}
#' \item{th}{p by p by nlam array of estimated interaction coefficients}
#' \item{obj}{nlam values of objective function, one per lambda value}
#' \item{lamlist}{Vector of values of lambda used}
#' \item{mx}{p-vector of column means of x}
#' \item{sx}{p-vector of column standard deviations of x}
#' \item{my}{mean of y}
#' \item{mzz}{column means of feature product matrix}
#' \item{szz}{column standard deviations of feature product matrix}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{predict.hierNet}, \link{hierNet.cv}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic.path(x,y)
#' print(fit)
#'
#' @export hierNet.logistic.path
hierNet.logistic.path <- function (x, y, lamlist=NULL, delta=1e-8, minlam=NULL, maxlam=NULL, flmin=.01, nlam=20,
diagonal=TRUE, strong=FALSE, aa=NULL,
zz=NULL, stand.main=TRUE, stand.int=FALSE,
rho=nrow(x), niter=100, sym.eps=1e-3,# ADMM params
step=1, maxiter=2000, backtrack=0.2, tol=1e-5, # GG params
trace=0) {
this.call=match.call()
stopifnot(y %in% c(0, 1))
x <- scale(x, center=TRUE, scale=stand.main)
mx <- attr(x, "scaled:center")
sx <- attr(x, "scaled:scale")
if (is.null(maxlam)) {
if (!is.null(minlam)) stop("Cannot have maxlam=NULL if minlam is non-null.")
maxlam <- max(abs(t(x) %*% y))
minlam <- maxlam * flmin
}
if (is.null(minlam)) minlam <- maxlam * flmin
if (is.null(lamlist))
lamlist <- exp(seq(log(maxlam), log(minlam), length=nlam))
nlam <- length(lamlist)
if (is.null(zz))
zz <- compute.interactions.c(x, diagonal=diagonal)
else stopifnot(is.matrix(zz))
zz <- scale(zz, center=TRUE, scale=stand.int)
mzz <- attr(zz, "scaled:center")
szz <- attr(zz, "scaled:scale")
zz <- as.numeric(zz)
p <- ncol(x)
cp2 <- choose(p, 2)
b0 <- rep(NA, nlam)
bp <- bn <- matrix(NA, nrow=p, ncol=nlam)
th <- array(NA, c(p, p, nlam))
obj <- rep(NA, nlam)
aa <- NULL
for (i in seq(nlam)) {
cat(c("i,lam=", i, round(lamlist[i],2)), fill=TRUE)
aa <- hierNet.logistic(x, y, lam=lamlist[i], delta=delta, diagonal=diagonal, strong=strong,
aa=aa, zz=zz, stand.main=FALSE, stand.int=FALSE,
rho=rho, niter=niter, sym.eps=sym.eps,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
trace=trace)
b0[i] <- aa$b0
bp[, i] <- aa$bp
bn[, i] <- aa$bn
th[, , i] <- aa$th
obj[i] <- aa$obj
}
dimnames(bp) <- dimnames(bn) <- list(as.character(1:p), NULL)
dimnames(th) <- list(as.character(1:p), as.character(1:p), NULL)
out <- list(b0=b0, bp=bp, bn=bn, th=th, obj=obj, lamlist=lamlist, delta=delta,
mx=mx, my=0, sx=sx, mzz=mzz, szz=szz,
type="logistic", diagonal=diagonal, strong=strong,
step=step, maxiter=maxiter, backtrack=backtrack, tol=tol,
call=this.call)
if (strong) {
# ADMM parameters:
out$rho <- aa$rho
out$niter <- niter
out$sym.eps <- sym.eps
}
class(out) <- "hierNet.path"
out
}
balanced.folds <- function(y, nfolds=min(min(table(y)), 10)) {
totals <- table(y)
fmax <- max(totals)
nfolds <- min(nfolds, fmax)
# makes no sense to have more folds than the max class size
folds <- as.list(seq(nfolds))
yids <- split(seq(y), y)
# nice we to get the ids in a list, split by class
###Make a big matrix, with enough rows to get in all the folds per class
bigmat <- matrix(NA, ceiling(fmax/nfolds) * nfolds, length(totals))
for(i in seq(totals)) {
bigmat[seq(totals[i]), i] <- sample(yids[[i]])
}
smallmat <- matrix(bigmat, nrow = nfolds) # reshape the matrix
### Now do a clever sort to mix up the NAs
smallmat <- permute.rows(t(smallmat)) ### Now a clever unlisting
# the "clever" unlist doesn't work when there are no N
# apply(smallmat, 2, function(x)
# x[!is.na(x)])
res <-vector("list", nfolds)
for(j in 1:nfolds) {
jj <- !is.na(smallmat[, j])
res[[j]] <- smallmat[jj, j]
}
return(res)
}
permute.rows <-function(x) {
dd <- dim(x)
n <- dd[1]
p <- dd[2]
mm <- runif(length(x)) + rep(seq(n) * 10, rep(p, n))
matrix(t(x)[order(mm)], n, p, byrow = TRUE)
}
#' Cross-validation function for hierNet
#'
#' Uses cross-validation to estimate the regularization parameter for hierNet
#'
#'
#' @param fit Object returned from call to hierNet.path or
#' hierNet.logistic.path. All parameter settings will be taken from this
#' object.
#' @param x A matrix of predictors, where the rows are the samples and the
#' columns are the predictors
#' @param y A vector of observations, where length(y) equals nrow(x)
#' @param nfolds Number of cross-validation folds
#' @param folds (Optional) user-supplied cross-validation folds. If provided,
#' nfolds is ignored.
#' @param trace Verbose output? 0=no, 1=yes
#' @return
#' \item{lamlist}{Vector of lambda values tried}
#' \item{cv.err}{Estimate of cross-validation error}
#' \item{cv.se}{Estimated standard error of cross-validation estimate}
#' \item{lamhat}{lambda value minimizing cv.err}
#' \item{lamhat.1se}{largest lambda value with cv.err less than or equal to
#' min(cv.err)+ SE }
#' \item{folds}{Indices of folds used in cross-validation}
#' \item{yhat}{n by nlam matrix of predicted values. Here, ith prediction is
#' based on training on all folds that do not include the ith data point.}
#' \item{nonzero}{Vector giving number of non-zero coefficients for each lambda
#' value}
#' \item{call}{The call to hierNet.cv}
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet},\link{hierNet.path},
#' \link{hierNet.logistic},\link{hierNet.logistic.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' fit=hierNet.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#' plot(fitcv)
#'
#'
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' y=1*(y>0)
#' fit=hierNet.logistic.path(x,y)
#' fitcv=hierNet.cv(fit,x,y)
#' print(fitcv)
#' plot(fitcv)
#'
#' @export hierNet.cv
hierNet.cv <- function(fit, x, y, nfolds=10, folds=NULL, trace=0) {
this.call <- match.call()
stopifnot(class(fit) == "hierNet.path")
if(fit$type=="gaussian"){errfun=function(y,yhat){(y-yhat)^2}}
if(fit$type=="logistic"){errfun=function(y,yhat){1*(y!=yhat)}}
n <- length(y)
if(is.null(folds)) {
folds <- split(sample(1:n), rep(1:nfolds, length = n))
}
else {
stopifnot(class(folds)=="list")
nfolds <- length(folds)
}
lamlist=fit$lamlist
# get whether fit was standardized based on fit$sx and fit$szz...
if (is.null(fit$mx)) stop("hierNet object was not centered. hierNet.cv has not been written for this (unusual) case.")
stand.main <- !is.null(fit$sx)
stand.int <- !is.null(fit$szz)
n.lamlist <- length(lamlist) ### Set up the data structures
size <- double(n.lamlist)
err2=matrix(NA,nrow=nfolds,ncol=length(lamlist))
for(ii in 1:nfolds) {
cat("Fold", ii, ":")
if(fit$type=="gaussian"){
a <- hierNet.path(x[-folds[[ii]],],y=y[-folds[[ii]]],
lamlist=lamlist, delta=fit$delta, diagonal=fit$diagonal, strong=fit$strong, trace=trace,
stand.main=stand.main, stand.int=stand.int,
rho=fit$rho, niter=fit$niter, sym.eps=fit$sym.eps, # ADMM parameters (which will be NULL if strong=F)
step=fit$step, maxiter=fit$maxiter, backtrack=fit$backtrack, tol=fit$tol) # GG descent params
yhatt=predict.hierNet(a,newx=x[folds[[ii]],])
}
if(fit$type=="logistic"){
a <- hierNet.logistic.path(x[-folds[[ii]],],y=y[-folds[[ii]]],
lamlist=lamlist, delta=fit$delta, diagonal=fit$diagonal, strong=fit$strong,
trace=trace, stand.main=stand.main, stand.int=stand.int,
rho=fit$rho, niter=fit$niter, sym.eps=fit$sym.eps, # ADMM parameters (which will be NULL if strong=F)
step=fit$step, maxiter=fit$maxiter, backtrack=fit$backtrack, tol=fit$tol) # GG descent params
yhatt=predict.hierNet.logistic(a,newx=x[folds[[ii]],])$yhat
}
temp=matrix(y[folds[[ii]]],nrow=length(folds[[ii]]),ncol=n.lamlist)
err2[ii,]=colMeans(errfun(yhatt,temp))
cat("\n")
}
errm=colMeans(err2)
errse=sqrt(apply(err2,2,var)/nfolds)
o=which.min(errm)
lamhat=lamlist[o]
oo=errm<= errm[o]+errse[o]
lamhat.1se=lamlist[oo & lamlist>=lamhat][1]
nonzero=colSums(fit$bp-fit$bn!=0) + apply(fit$th!=0, 3, function(a) sum(diag(a)) + sum((a+t(a)!=0)[upper.tri(a)]))
obj <- list(lamlist=lamlist, cv.err=errm,cv.se=errse,lamhat=lamhat, lamhat.1se=lamhat.1se,
nonzero=nonzero, folds=folds,
call = this.call)
class(obj) <- "hierNet.cv"
obj
}
plot.hierNet.cv <- function(x, ...) {
par(mar = c(5, 5, 5, 1))
yrang=range(c(x$cv.err-x$cv.se,x$cv.err+x$cv.se))
plot(log(x$lamlist), x$cv.err, xlab="log(lambda)",
ylab = "Cross-validation Error", type="n",ylim=yrang)
axis(3, at = log(x$lamlist), labels = paste(x$nonzero), srt = 90, adj = 0)
mtext("Number of features", 3, 4, cex = 1.2)
axis(2, at = c(0, 0.2, 0.4, 0.6, 0.8))
error.bars(log(x$lamlist), x$cv.err - x$cv.se, x$cv.err + x$cv.se, width = 0.01, col = "darkgrey")
points(log(x$lamlist), x$cv.err, col=2, pch=19)
abline(v=log(x$lamhat), lty=3)
abline(v=log(x$lamhat.1se), lty=3)
invisible()
}
error.bars <-function(x, upper, lower, width = 0.02, ...) {
xlim <- range(x)
barw <- diff(xlim) * width
segments(x, upper, x, lower, ...)
segments(x - barw, upper, x + barw, upper, ...)
segments(x - barw, lower, x + barw, lower, ...)
range(upper, lower)
}
#' Variable importance for hierNet.
#'
#' (This is an experimental function.) Calculates a measure of the importance
#' of each variable.
#'
#'
#' @param fit The results of a call to the "hierNet"
#' @param x The training set feature matrix used in call produced "fit"
#' @param y The training set response vector used in call produced "fit"
#' @param ... additional arguments (not currently used)
#' @return Table of variable importance.
#' @author Jacob Bien and Robert Tibshirani
#' @seealso \link{hierNet}, \link{hierNet.path}
#' @references Bien, J., Taylor, J., Tibshirani, R., (2013) "A Lasso for
#' Hierarchical Interactions." Annals of Statistics. 41(3). 1111-1141.
#' @examples
#'
#' set.seed(12)
#' x=matrix(rnorm(100*10),ncol=10)
#' x=scale(x,TRUE,TRUE)
#' y=x[,1]+2*x[,2]+ x[,1]*x[,2]+3*rnorm(100)
#' newx=matrix(rnorm(100*10),ncol=10)
#' fit=hierNet(x,y,lam=50)
#' yhat=predict(fit,newx)
#'
#' fit=hierNet.path(x,y)
#' yhat=predict(fit,newx)
#'
hierNet.varimp <- function(fit,x,y, ...) {
# NOTE: uses 0.5 cutoff for logistic case
lam=fit$lam
if(fit$type=="gaussian"){errfun=function(y,yhat){(y-yhat)^2}}
if(fit$type=="logistic"){
errfun=function(y,yhat){
term1=y*log(yhat);term1[yhat==0]=0
term2=(1-y)*log(1-yhat);term2[yhat==1]=0
val=-sum(term1+term2)
return(val)
}}
yhat=predict(fit,x)
rss=sum(errfun(y,yhat))
varsum=fit$bp-fit$bn+rowSums(abs(fit$th))
oo=which(abs(varsum)>1e-6)
imp=rss2=rep(NA,ncol(x))
for(j in oo){
cat(j)
fit0=fit;fit0$bp=fit$bp[-j];fit0$bn=fit$bn[-j];fit0$th=fit$th[-j,-j]
if(fit$type=="gaussian"){ fit2=hierNet(x[,-j],y,lam,delta=fit$delta,diagonal=fit$diagonal,aa=fit0)}
if(fit$type=="logistic"){ fit2=hierNet.logistic(x[,-j],y,lam,delta=fit$delta,diagonal=fit$diagonal,aa=fit0)}
yhat2=predict(fit2,x[,-j])
rss2[j]=sum(errfun(y,yhat2))
imp[j]=(rss2[j]-rss)/rss2[j]
}
imp[-oo]=0
res=cbind(1:ncol(x),round(imp,3))
ooo=order(-imp)
dimnames(res)=list(NULL,c("Predictor","Importance"))
cat("",fill=T)
return(res[ooo,])
}
Try the hierNet package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.