R/wirrls.R

In plsgenomics: PLS Analyses for Genomics

### wirrls.R  (2006-01)
###
###    Local Weighted Iteratively Reweighted Ridge Least Squares (IRRLS)
###                      for binary 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

wirrls <- function(Y, Z, Lambda=0, NbrIterMax=15, Threshold=10^(-12), WKernel) {
     ##    IN     
     ##########
     ##  Z : Design matrix
     ##          matrix n x (p+1)
     ##  Y : vector n 
     ##      response variable {0,1}-valued vector
     ##  Lambda :  coefficient of the ridge penalty
     ##          real
     ##  NbrIterMax : Maximal number of iterations
     ##          integer
     ##  Threshold : Used for the stopping rule.
     ##          real
     ##  WKernel : Kernel Weigth
     ##          matrix n x n 
     
     ##    OUT     
     ############
     ##  out :  structure that contains the fields
     ##      Gamma : vector of the regression coefficients w.r.t the design
     ##      matrix Z=[1 X].
     ##      Cvg : Cvg=1 if the algorithm has converged otherwise 0.
     
     dZ2 <- dim(Z)[2]
     n <- dim(Z)[1]
     p <- dZ2-1
     R <- matrix(0,dZ2,dZ2) 
     R[2:dZ2,2:dZ2] <- diag(1,nrow=p) 
     
     ##  THE NR ALGORITHM
     ##################################
     ##  1. Initialize the parameter
     
     c <- log(3) 
     Eta <- c*Y-c*(1-Y)                
     mu <- 1/(1+exp(-Eta))             
     DiagWNR <- mu*(1-mu)              
     WT <- diag(c(DiagWNR))%*%WKernel          
     H <- t(Z)%*%WT%*%Z+Lambda*R                
     trysolve <- try(solve(H, t(Z) %*% WKernel %*% ( diag(c(DiagWNR))%*%Eta+(Y-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   
          #Udapte mu
          mu <- 1/(1+exp(-Eta))   
          #Udapte Weight   
          DiagWNR <- mu*(1-mu)  
          WT <- diag(c(DiagWNR))%*%WKernel
          #Udapte Gradient
          Gradient <- t(Z)%*%WKernel%*%(Y-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((Eta>0)+0==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))
     return(list)
}