# R/rirls.R In lsplsGlm: Classification using LS-PLS for Logistic Regression

#### Documented in fit.lspcr.glm

```rirls <- function(Y,D,X,lambda=0,penalized=TRUE,nbrIterMax=15,threshold=10^(-12)) {
##    IN
##########
##  D : Design matrix
##          matrix n x q
##  X : Spectral design matrix
##          matrix n x p
##  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
##  penalized : if TRUE the parameter associated with D are ridge penalized

##    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.
##		W : matrix nxn weight at convergence
##		pseudoresponse : matrix nx1 pseudoresponse at convergence

p <- dim(X)[2]
q <- dim(D)[2]
n <- dim(X)[1]
Z <- cbind(rep(1,n),D,X)
dZ2 <- dim(Z)[2]

R<-matrix(0,dZ2,dZ2)
if(is.null(D)==TRUE)
{q<-0}
if (penalized==TRUE)
{R[2:dZ2,2:dZ2]<-diag(1,nrow=(p+q))}

R[(q+2):dZ2,(q+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))
H <- t(Z)%*%WT%*%Z+lambda*R
trysolve<-try(solve(H,t(Z)%*%(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
Eta[Eta>10]=10
Eta[Eta< -10]=-10
#Udapte mu
mu<-1/(1+exp(-Eta))
#Udapte Weight
DiagWNR <- mu*(1-mu)
WT <- diag(c(DiagWNR))
#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
StopRule <- max(Cvg,as.numeric(NbrIter>=nbrIterMax),Separation,Illcond)
}

##  CONCLUSION
###############

list <- list(coefficients=Gamma,cvg=max(Cvg,Separation),W=WT,pseudoresponse=(Eta+(Y-mu)/mu/(1-mu)))
return(list)
}

fit.lspcr.glm<-function(Y,X,D,ncomp){

k=ncomp
res <- prcomp(X,center=FALSE,scale=FALSE)
Xred <- res\$x[,1:k]
V <- res\$rotation[,1:k]
if (ncomp==1)
{
Xred <- matrix(Xred,ncol=1)
V <- matrix(V,ncol=1) ###  XLc%*%V-XLred &&&& t(V)%*%V=Id
}
Xpcr = cbind(D,Xred)
res=glm(Y~Xpcr,family = binomial)
reponse<-list(coefficients=res\$coefficients,cvg=res\$converged,projection=V)
}
```

## Try the lsplsGlm package in your browser

Any scripts or data that you put into this service are public.

lsplsGlm documentation built on July 27, 2017, 5:01 p.m.