R/selectVariables.R

   1 #' selectVaribles
   2 #' It is a function which construct, for a given lambda, the sets of
   3 #' relevant variables and irrelevant variables.
   4 #'
   5 #' @param phiInit an initial estimator for phi (size: p*m*k)
   6 #' @param rhoInit an initial estimator for rho (size: m*m*k)
   7 #' @param piInit    an initial estimator for pi (size : k)
   8 #' @param gamInit an initial estimator for gamma
   9 #' @param mini      minimum number of iterations in EM algorithm
  10 #' @param maxi      maximum number of iterations in EM algorithm
  11 #' @param gamma  power in the penalty
  12 #' @param glambda grid of regularization parameters
  13 #' @param X          matrix of regressors
  14 #' @param Y          matrix of responses
  15 #' @param thres  threshold to consider a coefficient to be equal to 0
  16 #' @param tau        threshold to say that EM algorithm has converged
  17 #'
  18 #' @return TODO
  19 #'
  20 #' @examples TODO
  21 #'
  22 #' @export
  23 selectVariables <- function(phiInit,rhoInit,piInit,gamInit,
  24     mini,maxi,gamma,glambda,X,Y,thres,tau)
  25 {
  26     dimphi <- dim(phiInit)
  27     p <- dimPhi[1]
  28     m <- dimPhi[2]
  29     k <- dimPhi[3]
  30     L <- length(glambda);
  31     A1 <- array(0, dim <- c(p,m+1,L))
  32     A2 <- array(0, dim <- c(p,m+1,L))
  33     Rho <- array(0, dim <- c(m,m,k,L))
  34     Pi <- array(0, dim <- c(k,L));
  35
  36     # For every lambda in gridLambda, comutation of the coefficients
  37     for (lambdaIndex in c(1:L))
  38     {
  39         Res <- EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,
  40             gamma,glambda[lambdaIndex],X,Y,tau);
  41         phi <- Res$phi
  42         rho <- Res$rho
  43         pi <- Res$pi
  44
  45         # If a coefficient is larger than the threshold, we keep it
  46         selectedVariables <- array(0, dim = c(p,m))
  47         discardedVariables <- array(0, dim = c(p,m))
  48         atLeastOneSelectedVariable <- false
  49         for (j in c(1:p))
  50         {
  51             cpt <- 1
  52             cpt2 <-1
  53             for (mm in c(1:m))
  54             {
  55                 if (max(abs(phi[j,mm,])) > thres)
  56                 {
  57                     selectedVariables[j,cpt] <- mm
  58                     cpt <- cpt+1
  59                     atLeastOneSelectedVariable <- true
  60                 } else
  61                 {
  62                     discardedVariables[j,cpt2] <- mm
  63                     cpt2 <- cpt2+1
  64                 }
  65             }
  66         }
  67
  68         # If no coefficients have been selected, we provide the zero matrix
  69         # We delete zero coefficients: vec = indices of zero values
  70         if (atLeastOneSelectedVariable)
  71         {
  72             vec <- c()
  73             for (j in c(1:p))
  74             {
  75                 if (selectedVariables(j,1) != 0)
  76                     vec <- c(vec,j)
  77                 # Else ( NOTE: [auder] else ?! TODO: explain? )
  78                 # we provide the indices of relevant coefficients
  79                 A1[,1,lambdaIndex] <- c(vec,rep(0,p-length(vec)))
  80                 A1[1:length(vec),2:(m+1),lambdaIndex] <- selectedVariables[vec,]
  81                 A2[,1,lambdaIndex] <- 1:p
  82                 A2[,2:(m+1),lambdaIndex] <- discardedVariables
  83                 Rho[,,,lambdaIndex] <- rho
  84                 Pi[,lambdaIndex] <- pi
  85             }
  86         }
  87     }
  88
  89     return(res = list(A1 = A1, A2 = A2 , Rho = Rho, Pi = Pi))
  90 }