R/mwirrls.R

In plsgenomics: PLS Analyses for Genomics

### mwirrls.R  (2006-01)
###
###    Local Weighted Iteratively Reweighted Ridge Least Squares (IRRLS)
###                      for categorical data
###
### Copyright 2006-01 Sophie Lambert-Lacroix and Julie Peyre
###
###
### This file is part of the `plsgenomics' library for R and related languages.
### It is made available under the terms of the GNU General Public
### License, version 2, or at your option, any later version,
### incorporated herein by reference.
### 
### This program is distributed in the hope that it will be
### useful, but WITHOUT ANY WARRANTY; without even the implied
### warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
### PURPOSE.  See the GNU General Public License for more
### details.
### 
### You should have received a copy of the GNU General Public
### License along with this program; if not, write to the Free
### Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
### MA 02111-1307, USA

mwirrls <- function(Y, Z, Lambda=0, NbrIterMax=15, Threshold=10^(-12), WKernel) {
     ##    IN     
     ##########
     ##      c=NbrClass-1
     ##  Z : Design matrix (cn) x c(p+1)
     ##          block data matrix
     ##  Y : Response Vector
     ##          vector ntrain 
     ##  Lambda :  coefficient of the ridge penalty
     ##          real
     ##  NbrIterMax : Maximal number of iterations
     ##          positive integer
     ##  Threshold : Used for the stopping rule.
     ##          real
     ##  WKernel : Kernel Weigth
     ##          matrix cn x cn 
     
     ##    OUT     
     ############
     ##      out :  structure that contains the fields
     ##      Gamma : vector of the regression coefficients w.r.t the design
     ##      matrix.
     ##      Cvg : Cvg=1 if the algorithm has converged otherwise 0.
     ##      Ybloc : block Response Vector ntrain*c
     
     c <- max(Y)
     dZ2 <- dim(Z)[2]
     n <- length(Y)
     p <- dZ2/c-1
     
     Ybloc <- rep(0,(n*c))
     for (cc in 1:c) {
          ff <- which(Y==cc)  
          Ybloc[(ff-1)*c+cc]=rep(1,length(ff))
     }
     
     R <- matrix(0,dZ2,dZ2) 
     diag(R) <- rep(1,dZ2)
     grid <- seq(from=1, to=dZ2, by= (p+1))
     R[grid,grid] <- 0 
     
     
     
     ##  THE NR ALGORITHM
     ##################################
     ##  1. Initialize the parameter
     
     mu <- (1+2*Ybloc)/(c+3)
     Eta <- rep(0,length(mu))
     WNR <- matrix(0,length(mu),length(mu))
     
     for (kk in 1:n) {
          ss <- 1-sum(mu[c*(kk-1)+(1:c)])
          Blocmu <- mu[c*(kk-1)+(1:c)]
          BlocW <- -Blocmu%*%t(Blocmu)
          BlocW <- BlocW+diag(Blocmu)
          WNR[c*(kk-1)+(1:c),c*(kk-1)+(1:c)] <- BlocW
          Eta[c*(kk-1)+(1:c)] <- log(Blocmu)-log(ss)
     }
     
     WT <- WKernel%*%WNR           
     H <- t(Z)%*%WT%*%Z+Lambda*R    
     trysolve<-try(solve(H,t(Z)%*%WKernel%*%(WNR%*%Eta+(Ybloc-mu))), silent = TRUE)
     if (sum(is.nan(trysolve))>0) {
          trysolve <- NULL
     }
     
     if (is.matrix(trysolve)==FALSE) {
          Gamma <- rep(1,dZ2) 
          StopRule <- 1   
          NbrIter <- NbrIterMax+1
          Illcond <- 1
          Separation <- 0
          Cvg <- 0
     }
     
     if (is.matrix(trysolve)==TRUE) {
          StopRule <- 0 
          NbrIter <- 0
          Illcond <- 0
          Separation <- 0
     }              
     
     ##  2. Newton-Raphson loop
     
     while (StopRule==0) {
          
          #Increment the iterations number 
          NbrIter <- NbrIter + 1
          
          #Udapte Gamma
          if (NbrIter==1) {
               Gamma <- trysolve
          }
          if (NbrIter>1) {
               Gamma <- Gamma+trysolve
          }
          
          #Udapte Eta            
          Eta <- Z%*%Gamma  
          Eta[Eta>700] <- rep(700,sum(Eta>700)) 
          
          #Udapte mu and WNR
          for (kk in 1:n) {
               mu[c*(kk-1)+(1:c)] <- exp(Eta[c*(kk-1)+(1:c)])/(1+sum(exp(Eta[c*(kk-1)+(1:c)])))
               Blocmu <- mu[c*(kk-1)+(1:c)]
               BlocW <- -Blocmu%*%t(Blocmu)
               BlocW <- BlocW+diag(Blocmu)
               WNR[c*(kk-1)+(1:c),c*(kk-1)+(1:c)] <- BlocW
          }
          
          #Udapte total Weight   
          WT <- WNR%*%WKernel
          
          #Udapte Gradient
          Gradient <- t(Z)%*%WKernel%*%(Ybloc-mu)-Lambda*R%*%Gamma
          
          #Udapte H           
          H <- t(Z)%*%WT%*%Z+Lambda*R
          trysolve <- try(solve(H,Gradient), silent = TRUE)
          if (sum(is.nan(trysolve))>0) {
               trysolve <- NULL
          }
          
          #Compute the StopRule
          #on the (quasi)-separation detection (for Lambda=0)
          if (Lambda==0) {
               Separation <- as.numeric(sum((apply(cbind(rep(0,n),matrix(Eta,nrow=n,byrow=TRUE)),1,which.max)-1)==Y)==n)
          }
          
          #on the conditioning of matrix
          Illcond <- as.numeric(is.matrix(trysolve)==FALSE)
          
          #on the convergence
          Cvg <- as.numeric(sqrt(sum(abs(Gradient)^2))<=Threshold)
          StopRule <- max(Cvg,as.numeric(NbrIter>=NbrIterMax), Separation, Illcond)
          
     }
     
     ##  CONCLUSION
     ###############        
     
     list <- list(Coefficients=Gamma, Cvg=max(Cvg,Separation), Ybloc=Ybloc)
     return(list)
}