R_MM_backup/NIPALS_PLS_component.R

In kimsche/MetaboMate: All you need for NMR-based metabolic profiling with R!

#' Calculating a single PLS component
#' @param X Input matrix with rows and columns representing observations and variables
#' @param Y Dependend variable, in form of dummy matrix (multi-levels allowed) or numeric column vector
#' @return Returned is a list with the following entries:
#' \item{X.res}{Residual X matrix.}
#' \item{Y.res}{Residual Y matrix.}
#' \item{scores}{PLS component scores.}
#' \item{loadings}{PLS component loadings.}
#' \item{weights}{PLS component weights.}
#' \item{betas}{PLS X coefficients.}
#' \item{Q.pc}{PLS Y coefficient.}
#' @author Torben Kimhofer \email{tkimhofer@@gmail.com}
#' @noRd

NIPALS_PLS_component<-function(X, Y){

  # X=T*P +E
  # Y=U*Q +F*

  dd<-1
  count<-1

  # initialise scores (t_h)
  u_h <- cbind(Y[,1]) # u_h start is Y
  while(dd>1e-10){

    # X block: calc weights, normalise and cacultate scores (t_h)
    w_h <- t(X) %*% u_h / drop(crossprod(u_h))
    w_h <- w_h / sqrt(sum(w_h^2)) # scale to length 1

    # calc X scores (t_h)
    t_h <- X %*% w_h / drop(crossprod(w_h)) # normalisation not needed (is one)

    # check convergence of t with last iteration (dd crit in while loop)
    if(count>1){
      dd<-sum((t_h[,1]-t_old[,1])^2);
    }
    t_old<-t_h

    ## Y block, calc weights, normalise and calc Y scores
    ## steps can be omitted for 2 class Y (simply by setting q_h=1)
    q_h<-t(Y) %*% t_h /  drop(crossprod(t_h))
    q_h<-q_h/sqrt(sum(q_h^2)) # normalisation

    u_h<-(Y) %*% (q_h) / drop(crossprod((q_h)))
    count<-count+1
  }

    # calculate the X loadings and rescale scores and weights accordingly
    p_h <- (t(X) %*% t_h) / drop(crossprod(t_h))
    p_h <- p_h / sqrt(sum(q_h^2))

    # calc y loadings
    yl<-t(Y) %*% u_h/drop(crossprod(u_h))

    # calc b for residuals Y
    # calculate beta (regression coefficient) via inverse insted of subtracting q_h directly
    b<-(t(u_h) %*% t_h) /drop(crossprod(t_h))

    # # calc b for residuals Y
    # # calculate beta (regression coefficient) via inverse insted of subtracting q_h directly
    # b= solve(t(t_h) %*% t_h) %*% t(t_h) %*% cbind(u_h)
    # Y_res = Y - t( b[1,1] *t(t_h)) %*% (q_h)

  # calc residual matrice X and Y
  X_res <- X - (t_h %*% t(p_h))
  Y_res <- Y - (drop(b) *t_h) %*% t(q_h)


  # define slots for PLS_Torben object
  res<-list(X.res=X_res, Y.res=Y_res, scores=t_h, loadings=t(p_h), weights=t(w_h), betas=as.numeric(b), Qpc=q_h)
  #
  # track = setClass('PLS_Torben', slots=c(Xres='matrix', Yres='matrix', scores='matrix', loadings='matrix', weights='matrix', betas='numeric', Qpc='matrix')) #, Xcenter='numeric', Xscale='numeric', Ycenter='numeric', Yscale='numeric'
  #
  # # create slots for PLS_Torben objects
  # mod_pls=track(Xres=X_res, Yres=Y_res, scores=t_h, loadings=t(p_h), weights=t(w_h), betas=as.numeric(b), Qpc=q_h) #, Xcenter=meanX, Xscale=sdX, Ycenter=meanY, Yscale=sdY
  #

  return(res)
}