R/HCgglasso.R

In HCgglasso: Group-Lasso with Hierarchical Clustering

#' Run hierarchical clustering following by a group-lasso on all the different partition.
#'
#' @author Quentin Grimonprez
#' @param X matrix of size n*p
#' @param y vector of size n
#' @param hc output of \code{\link{hclust}} function. If not provided, \code{\link{hclust}} is run with ward.D2 method
#' @param lambda lambda values for group lasso. If not provided, the function generates its own values of lambda
#' @param weightLevel a vector of size p for each level of the hierarchy. A zero indicates that the level will be ignored. If not provided, use 1/(height between 2 successive levels)
#' @param weightSizeGroup a vector of size 2*p-1 containing the weight for each group. Default is the square root of the size of each group
#' @param intercept should an intercept be included in the model ?
#' @param verbose print some information
#' @param ... Others parameters for \code{\link{gglasso}} function
#'
#' @return a HCgglasso object containing :
#' \describe{
#'   \item{lambda}{lambda values}
#'   \item{b0}{intercept values for \code{lambda}}
#'   \item{beta}{A list containing the values of estimated coefficients for each values of \code{lambda}}
#'   \item{var}{A list containing the index of selected variables for each values of \code{lambda}}
#'   \item{group}{A list containing the values index of selected groups for each values of \code{lambda}}
#'   \item{nVar}{A vector containing the number of non zero coefficients for each values of \code{lambda}}
#'   \item{nGroup}{A vector containing the number of non zero groups for each values of \code{lambda}}
#'   \item{structure}{A list containing 3 vectors. var : all variables used. group : associated groups. 
#'   weight : weight associated with the different groups.
#'   level : for each group, the corresponding level of the hierarchy where it appears and disappears. 3 indicates the level with a partition of 3 groups.}
#'   \item{time}{computation time}
#'   \item{dim}{dimension of \code{X}}
#'   \item{hc}{Output of hierarchical clustering}
#'   \item{call}{Code executed by user}
#' } 
#'
#'
#' @examples 
#' set.seed(42)
#' X = simuBlockGaussian(50,12,5,0.7)
#' y = drop(X[,c(2,7,12)]%*%c(2,2,-2)+rnorm(50,0,0.5))
#' res = HCgglasso(X,y)
#' 
#' @seealso \link{cv.HCgglasso}, \link{stability.HCgglasso}, \link{listToMatrix}, \link{predict.HCgglasso}, \link{coef.HCgglasso}, \link{plot.cv.HCgglasso}
#' 
#' @export
HCgglasso <- function(X, y, hc = NULL, lambda = NULL, weightLevel = NULL, weightSizeGroup = NULL, intercept = TRUE, verbose = FALSE,...)
{
  #check parameters 
  .checkParameters(X, y, hc, lambda, weightLevel, weightSizeGroup, intercept, verbose)
  
  
  # define some usefull variables
  n = nrow(X)
  p = ncol(X)
  tcah = rep(NA,3)
  
  ######## hierarchical clustering
  #if no hc output provided, we make one
  if(is.null(hc))
  {
    if(verbose)
      cat("Computing hierarchical clustering...")
    
    t1 = proc.time()
    d = dist(t(X))
    hc = fastcluster::hclust(d, method = "ward.D2")
    t2 = proc.time()
    tcah = t2-t1
    if(verbose)
      cat("DONE in ",tcah[3],"s\n")
  }
  
  ######## compute weight, active variables and groups
  if(verbose)
    cat("Preliminary step...")
  t1 = proc.time()
  prelim = preliminaryStep(hc, weightLevel, weightSizeGroup)  


  #duplicate data
  Xb = X[,prelim$var]
  t2 = proc.time()
  if(verbose)
    cat("DONE in ",(t2-t1)[3],"s\n")
  
  ######## group lasso
  if(verbose)
    cat("Computing group-lasso...")
  t1 = proc.time()
  res = gglasso(Xb,y,prelim$group,pf=prelim$weight,lambda=lambda,intercept=intercept,...)
  t2 = proc.time()
  tgglasso = t2-t1
  if(verbose)
    cat("DONE in ",tgglasso[3],"s\n")
  
  ########  create output object
  res2 = list()
  res2$lambda = res$lambda
  non0 = apply(res$beta,2,FUN=function(x){which(x!=0)})
  res2$var = lapply(non0,FUN=function(x){prelim$var[x]})
  res2$nVar = sapply(res2$var,FUN=function(x){length(unique(x))})
  res2$group = lapply(non0,FUN=function(x){prelim$group[x]})
  res2$nGroup = sapply(res2$group,FUN=function(x){length(unique(x))})
  res2$beta = lapply(1:length(res$lambda),FUN=function(x){res$beta[non0[[x]],x]})
  res2$b0 = res$b0
  res2$structure = prelim
  res2$dim = dim(X)
  res2$hc = hc
  res2$time = c(tcah[3],tgglasso[3])
  res2$call = match.call()
  names(res2$time) = c("hclust","glasso")    
  class(res2) = "HCgglasso"
  
  
  return(res2)
}


#
# compute the mnimimum weight of each group
# 
# @param hc outup of hclust function
#
levelMinWeight <- function(hc, weightLevel = NULL)
{
  p=length(hc$order)
  
  # highest level at which cluster are seen for the last time
  # the p first are the single variables and the p-2 next are the cluster in the order of apparition
  lvSingle = sapply((-1):(-p),FUN=function(i){which(hc$merge==i)%%(p-1)})
  lvCluster = sapply(1:(p-2),FUN=function(i){which(hc$merge==i)%%(p-1)})
  lvCluster[lvCluster==0] = p-1
  lvCluster = c(lvCluster,p)
  
  # branch length. The first one is associated with the partition in 2 clusters
  if(is.null(weightLevel))
    weightLevel = c(0,sqrt(1/diff(hc$height)),0)
    
  # minimum weight of levels of each cluster
  minLevelWeight = rep(0, 2*p-1)
  # minimal weight of single variable
  minLevelWeight[1:p] = sapply(lvSingle,FUN=function(i){ind=(weightLevel[1:i]!=0);ifelse(sum(ind),min(weightLevel[1:i][ind]),0)})#If there is only 0, we return 0, else we return the min > 0
  # mininmal weight for groups of 2 and more variables
  minLevelWeight[(p+1):(2*p-1)] = sapply(1:length(lvCluster),FUN=function(i){ind=(weightLevel[(i+1):lvCluster[i]]!=0);ifelse(sum(ind),min(weightLevel[(i+1):lvCluster[i]][ind]),0)})
  
  return(minLevelWeight)
}

#
# @param hc output of hierarchical clustering
#
# @return A matrix with 2 rows, the first row contains the level at which appears each group during the hierarchical clustering
# the second row contains the last level where the group is present. The p first columns represent single variable, the other the cluster in the order
# they appear in the hierarchical clustering
#
levelGroupHC <- function(hc)
{
  # Number of variables in the HC
  p <- nrow(hc$merge) + 1
  
  # Output matrix
  startend <- matrix(nrow = 2, ncol = 2*p-1)
  # Level where first appeared each group (j = level containg j groups)
  startend[1,] = c(rep(p,p),(p-1):1)
  
  # Find the level where each group diassapear
  for(i in 1:nrow(hc$merge))
  {
    for(j in 1:2)
    {   
      # a negative number indicates a single variable
      if(hc$merge[i,j] < 0)
      {
        startend[2,abs(hc$merge[i,j])] = p - i + 1
      }
      else # a positive number indicates a cluster of 2 or more variables
      {
        startend[2,p + hc$merge[i,j]] = p - i + 1
      }
      
    }# end for col os hs$merge
  }# end for row of hc$merge
  
  # Last group containing all variables
  startend[2,ncol(startend)] = 1
  
  rownames(startend) <- c("start", "end")
  
  return(startend)
}


#
# preliminary step for HCgglasso. Compute weight, active variables and groups
#
preliminaryStep <- function(hc, weightLevel = NULL, weightSizeGroup = NULL)
{
  #find unique groups of the hclust output
  uni = uniqueGroupHclust(hc)
  
  ######## Compute weights
  #compute the minimal weight of partition
  weightLevelGroup = levelMinWeight(hc, weightLevel)
  
  # CORRECTION : If weight is infinite, we change in 0 and it will be ignored
  weightLevelGroup[which(is.infinite(weightLevelGroup))] = 0
  
  #weight for group size
  if(is.null(weightSizeGroup))
    weightSizeGroup = as.vector(sqrt(table(uni$indexGroup)))
  
  #weight for each group
  weight = weightSizeGroup*weightLevelGroup
  
  # new weight without ignored groups
  weightb = weight
  ignoredGroup = which(weight==0)#groups with 0 weights
  weightb = weightb[-ignoredGroup]#we delete zeros weight
  
  #level of hc associated to groups
  p = length(hc$order)
  lv = levelGroupHC(hc)
  lv = lv[,-ignoredGroup]
  
  ######## Create data for gglasso
  varToDelete = uni$indexGroup%in%ignoredGroup
  var = uni$varGroup[!varToDelete]
  group = uni$indexGroup[!varToDelete]
  
  #group must be consecutively numbered 1,2,3,...
  #need a correction when some groups have to be ignored
  if(length(ignoredGroup)>0)
  {
    difNumber = rep(0,length(group))
    for(i in 1:length(ignoredGroup))
    {
      ind = which(group>ignoredGroup[i])
      difNumber[ind] = difNumber[ind]-1
    }
    group = group + difNumber
  }
  
  return(list(group=group,var=var,weight=weightb,level=lv))
}

# check parameters of HCgglasso function
.checkParameters <- function(X, y, hc, lambda, weightLevel, weightSizeGroup, intercept, verbose)
{
  #check X
  if(!is.matrix(X)) 
    stop("X has to be a matrix.")
  if(any(is.na(X))) 
    stop("Missing values in X not allowed.")
  if(!is.numeric(X))
    stop("X has to be a matrix of real.")
  
  #check y
  if(!is.numeric(y))
    stop("y has to be a vector of real.")
  if(any(is.na(y))) 
    stop("Missing values in y not allowed.")
  
  #check if X and y are compatible
  if(nrow(X)!=length(drop(y)))
    stop("The length of y and the number of rows of X don't match.")
  
  #check hc
  if(!is.null(hc))
  {
    #check if hc is a hclust object
    if(class(hc)!="hclust")
      stop("hc must be an hclust object.")
    #check if hc and X are compatible
    if(length(hc$order)!=ncol(X))
      stop("hc is not a clustering of the p covariates of X.")
    
  }
  
  #check if lambda is a vector of positive real
  if(!is.null(lambda))
  {
    if(!is.numeric(lambda))
      stop("lambda must be a vector of positive real.")
    if(any(lambda < 0))
      stop("lambda must be a vector of positive real.")
  }
  
  #check if weightLevel is a vector of positive real
  if(!is.null(weightLevel))
  {
    if(!is.numeric(weightLevel))
      stop("weightLevel must be a vector of positive real.")
    if(length(weightLevel)!=ncol(X))
      stop("weightLevel must have the same length as the number of columns of matrix X.")
    if(any(weightLevel < 0))
      stop("weightLevel must be a vector of positive real.")
  }
  
  #check if weightSizeGroup is a vector of positive real
  if(!is.null(weightSizeGroup))
  {
    if(!is.numeric(weightSizeGroup))
      stop("weightSizeGroup must be a vector of real.")
    if(any(weightSizeGroup < 0))
      stop("weightSizeGroup must be a vector of positive real.")
  }
  
  #check if intercept is a boolean
  if(length(intercept)!=1)
    stop("intercept must be a boolean.")
  if(!is.logical(intercept))
    stop("intercept must be a boolean.")
  
  #check if verbose is a boolean
  if(length(verbose)!=1)
    stop("verbose must be a boolean.")
  if(!is.logical(verbose))
    stop("verbose must be a boolean.")
  
  invisible(return(NULL))
}