R/regression_functions.R

In WarrenDavidAnderson/dynamicNetworkID:

#' Y vector for regression
#'
#' This function establishes the Y vector for a specific set of m-values (Equation 11, Anderson et al., 2017).
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for 
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). Defaults to datHMF.
#' @param Y The Y data obtained with getXY(). Defaults to NULL.
#' @param numPar The number of parameters found by multiplying the number of annotations 
#' by the number of measured features (e.g., number of orggans * number of genes). 
#' numPar = ((ncol(datHMF)-2)^2)*(length(unique(datHMF[,2]))^2). Defaults to NULL.
#' @param mRange The m-value range. Defaults to NULL.
#' @return The Y vector for regression analysis
#' @export
#' @examples Yreg = get_Yreg(datHMF,XY$Y,numPar,mRange)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
get_Yreg = function(datHMF=NULL,Y=NULL,numPar=NULL,mRange=NULL){
  Yreg <- matrix( c(0), sqrt(numPar)*length(mRange), 1 )
  x = 1
  len = length(mRange)
  for (annots in 1:length(unique(datHMF[,2]))) {
    for (feats in 1:(ncol(datHMF)-2)) {
      rowsY <- ( length(mRange)*(x-1)+1 ) : ( length(mRange)*x )
      Yreg[rowsY] <- Y[feats,annots,mRange]
      x = x + 1
    } 
  }
  return(Yreg)
} # end get_Yreg



#' X matrix for regression
#'
#' This function establishes the X matrix for a specific set of m-values (Equation 11, Anderson et al., 2017).
#' @param datHMF The data matrix with scaled measurements measurements. Scaled center measurements are provided for 
#' each measurement/annotation (e.g., gene/organ) combination. This format is produced by scale_zeroOne(). Defaults to datHMF.
#' @param X The X data obtained with getXY(). Defaults to NULL.
#' @param numPar The number of parameters found by multiplying the number of annotations 
#' by the number of measured features (e.g., number of orggans * number of genes). 
#' numPar = ((ncol(datHMF)-2)^2)*(length(unique(datHMF[,2]))^2). Defaults to NULL.
#' @param mRange The m-value range. Defaults to NULL.
#' @return The Y vector for regression analysis
#' @export
#' @examples Xreg = get_Xreg(datHMF,XY$Y,numPar,mRange)
#' @references \url{http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005627}
get_Xreg = function(datHMF=NULL,X=NULL,numPar=NULL,mRange=NULL){
  Nr = length(unique(datHMF[,2])) # number of organs (annotations)
  Xreg <- matrix( c(0), sqrt(numPar)*length(mRange), numPar  )
  x = 1
  for (annots in 1:Nr) {
    for (feats in 1:(ncol(datHMF)-2)) {
      
      # all genes in all organs influence gene g in organ o:
      XX <- matrix(c(0),length(mRange),sqrt(numPar))
      indXX <- 1
      for (an in 1:length(unique(datHMF[,2]))) {
        for (ft in 1:(ncol(datHMF)-2)) {
          XX[,indXX] <- X[ft,an,mRange] 
          indXX <- indXX + 1
        } 
      } 
      
      # add data to the X regression matrix
      rowsX <- ( length(mRange)*(x-1)+1 ) : ( length(mRange)*x ) 
      colsX <- ( sqrt(numPar)*(x-1)+1 ) : ( sqrt(numPar)*x )
      Xreg[rowsX,colsX] <- XX
      x <- x + 1
      
    } # feats
  } # annots
  return(Xreg)
} # end get_Xreg


#' Get rate constants
#'
#' This function formats the rate constant data after the fir has been implemented using glmnet() from the glmnet library. 
#' Sometimes glmnet() produces less fits than requested based on the nlambda parameter. In such cases, get_rateConst() 
#' pads the parameter matrix with zeros such that the size of the matrix reflects the selection of nlambda.
#' @param fit The beta values from the fit object generated by glmnet().
#' @param nlambda The number of lambda values.
#' @return The rate constant matrix.
#' @examples rateConst = get_rateConst(fit$beta,nlambda)
#' @references \url{https://www.jstatsoft.org/article/view/v033i01}
get_rateConst = function(fit=NULL,nlambda=NULL){
  rateConst <- data.matrix(fit) 
  if(ncol(rateConst)<nlambda){
    del = nlambda - ncol(rateConst)
    addcol = matrix(c(0),nrow(rateConst),del)
    rateConst = cbind(rateConst,addcol)
  }
  return(rateConst)
} # end get_rateConst