R/LoBrA.R

In LoBrA: Generalized Spline Mixed Effect Models for Longitudinal Breath Data

#' LoBrA: A package for modeling longitudinal breath data
#'
#' The \code{LoBrA} package provides important data objects and functions to analyze longitudinal metabolomic (breath) data.
#' 
#' @section Introduction:
#' Novel metabolomic technologies paved the way for longitudinal analysis of exhaled air 
#' and online monitoring of fast progressing diseases. 
#' This package implements an automated analysis approach of longitudinal data from 
#' different omics technologies, such as ion mobility spectrometry of human exhaled air 
#' and demonstrates how including temporal signals increases the statistical power in 
#' biomarker identification. 
#' It can handel multiple irregular 4D time series data. More precisely, it can simultaniously 
#' handel the data of multiple experiements each observing multiple components. 
#' Therefore, it allows repeated measurements of a component, irregular sampling, and unequal 
#' temporal spacing of the time points.
#' 
#' @section LoBrA Analysis:
#' A typical LoBrA analysis is will comprise the following steps
#' 
#' 1. Background Screening: Using the function \code{\link{screening}} and \code{\link{selectComponents}} to select the Components that most likely do not originate from background noise.
#' 
#' 2. Model Selection: First, a set of spline models based on different number of splits and split positions are generated by the function \code{\link{lobraModelSelection}}. Subsequently, these models are evaluated using different quality criteria, i.e. 'AIC', 'BIC' and 'LogLik'. Finally, the most appropriate model is selected.
#' 
#' 3. Evaluation of the non-background components on the selected model, using the longitudinal 'Gouderman' linear mixed effect model in function \code{\link{modelGoudermanLongitudinal}}. 
#' 
#' @section Author(s):
#' Maintainer:
#' Anne-Christin Hauschild [Copyright holder]
#' 
#' Authors:
#' \itemize{
##'  \item{Sandrah P. Eckel}
##'  \item{Jan Baumbach}
##' }
#' 
#'
#' @docType package
#' @name LoBrA
#' 
NULL


#' An S4 class to represent a 'LoBrA' Data Object (LDO). It stores multiple time series data for muliple experiements and multiple Components. It allows repeated measurements of a component, irregular sampling, and unequal temporal spacing of the time points.
#'
#' @slot name	character Name of the 'LDO' object
#' @slot dataMatrices	list List of matrices of component measurement. It contains a measurement for each time point and each experiment.
#' @slot backgroundMatrices	list List of matrices of background measurements. It contains a measurement for each time point and each experiment.
#' @slot peaknames	character Character vector of Component names
#' @slot times	numeric Vector of times for each time point in the data
#' @slot ids	character Vector of identifiers for the experiments
#' @slot labels	factor Vector of class labels for each experiment
#' @importFrom methods setClass 
#' @export 
#' 
LDO <- methods::setClass("LDO", slots=list(name="character", dataMatrices="list", backgroundMatrices="list", peaknames="character", times="numeric", ids="character", labels="factor"))

#' An S4 class to represent the screening of metabolites in an 'LDO'.
#' 
#' @slot ldo	LDO LDO object the screening is based on.
#' @slot experimentIntercept	list List of experiment intercepts.
#' @slot experimentResiduals	list List of experiment Residuals.
#' @slot interceptPvalues	matrix Matrix of experiment intercept p-values.
#' @slot residualPvalues	matrix Matrix of experiment Residual p-values.
#' @slot selectedPeaks	matrix Matrix of logical values. Each entry indicates whether a specific component is significant according to a specific test.
#' @importFrom methods setClass 
#' @export 
#'
LDOscreening<- methods::setClass("LDOscreening", slots=list(ldo="LDO", experimentIntercept="list", experimentResiduals="list",  interceptPvalues="matrix", residualPvalues="matrix", selectedPeaks="matrix"))


#' An S4 class to represent a model selection result based on an 'LDO'.
#' @slot ldo	LDO 'LDO' object the model selection is based on.
#' @slot potentialBreaks	numeric Vector of numeric values that were considered as potential break points in the model selection.
#' @slot splinetype	character Type of spline used.
#' @slot qualityMeasure	character Quality measures used during the model selection ('AIC', 'BIC' or 'LogLik')
#' @slot modelList	list List of all models tested.
#' @slot quality	list List of quality matrices, one matrix for each quality measure used. Each matrix contains the quality for each spline tested for each component.
#' @slot breaks	list For each tested spline, this list contains a vector of breaks.
#' @importFrom methods setClass 
#' @export 
#'
LDOmodelselection<- methods::setClass("LDOmodelselection", slots=list(ldo="LDO", potentialBreaks="numeric", splinetype="character", qualityMeasure="character", modelList="list", quality="list", breaks="list"))



#' An S4 class to represent a 'Gouderman' LDO object, that was generated by the generalized gauderman algorithm.
#' @slot name	character Name of the new 'generalized-Gauderman' adjusted longitudinal data
#' @slot dataFrames	list List of 'generalized-Gauderman' modified data. One \code{data.frame} for each component.
#' @slot peaknames	character Vector of component names contained in this object.
#' @slot k	numeric Updated times for the breaks of the spline model.
#' @slot times	matrix Vector of updated time values.
#' @slot newTimeVars	character The names of the newly defined time variables of the generalized 'Gauderman' model.
#' @slot ids	character Vector of identifiers for the experiments
#' @slot labels	factor Vector of class labels for each experiment
#' @importFrom methods setClass 
#' @export 
#'
GaudermanLDO <- methods::setClass("GaudermanLDO", slots=list(name="character", dataFrames="list", peaknames="character", k="numeric", times="matrix", newTimeVars="character", ids="character", labels="factor"))

#' An S4 class to represent the result of the linear mixed effect modeling on a gauderman LDO.
#' @slot name	character Name of the new 'generalized-Gauderman' adjusted longitudinal model.
#' @slot gaudermanLDO	GaudermanLDO 'Generalized-Gauderman' adjusted longitudinal data object.
#' @slot models	list List of models generated for each component.
#' @slot labels factor Vector of class labels for each experiment
#' @slot pvalues	matrix Matrix of p-values for the intercept as well as all slops of the spline model for each component.
#' @slot correctedpvalues	matrix Matrix of corrected p-values for the intercept as well as all slops of the spline model for each component.
#' @slot modelparameter	matrix Model parameter for each component.
#' @importFrom methods setClass 
#' @export 
GaudermanModelEvaluation <- methods::setClass("GaudermanModelEvaluation", slots=list(name="character", gaudermanLDO="GaudermanLDO", models="list", labels="factor", pvalues="matrix", correctedpvalues="matrix", modelparameter="matrix"))



#' @title Transformation of a single longitudinal data matrix into 'LoBrA' Data Object.
#' 
#' @description Real signals and background noise originating from experimental settings or random events 
#' @param longData Matrix of longitudinal data containing all components
#' @param name name of the dataset
#' @param id name to identify the experiment id column
#' @param time name to identify the time column
#' @param type name to identify the type column
#' @param class name to identify the class column
#' @param bg indicates whether the data table contains background data
#' @importFrom methods new
#' @export
#' @return 'LoBrA' data object
#' @examples \dontrun{
#' } 
#'   data(LoBraExample)
#'   name="Longitudinal Test Dataset"
#'   ldo<-as.LOBdataset(longDataExample, name, bg=TRUE)
#'   
as.LOBdataset<-function(longData, name="", id="id", time="time", type="type", class="class", bg=FALSE){

  peaknames<-colnames(longData)
  peaknames<- peaknames[(peaknames!=id)&(peaknames!=time)&(peaknames!=type)&(peaknames!=class)]
  times=sort(unique(longData[,time]))
  ids=sort(unique(longData[,id]))
  labels=as.factor(sort(unique(longData[,class])))
  bgmatrix<-NA;
  datamatrix<-NA;
  if(bg){
    bgmatrix<-list();
    datamatrix<-list();
    bdata<-longData[longData[,type]=="b",];
    sdata<-longData[longData[,type]=="s",];
    # dim(bdata)
    # dim(sdata)
    # dim(longData)
    peak<-peaknames[1];
    for(peak in peaknames){
      newSampleMatrix <- data.frame(class=factor(sdata[,class], levels=labels), 
                                   id = factor(sdata[,id], levels=ids), 
                                   time=as.numeric(sdata[,time]), 
                                   value=as.numeric(sdata[,peak]))
      newbackgroundMatrix <- data.frame(class=factor(bdata[,class], levels=labels), 
                                   id = factor(bdata[,id], levels=ids), 
                                   time=as.numeric(bdata[,time]), 
                                   value=as.numeric(bdata[,peak]))
      bgmatrix[[peak]]<-newbackgroundMatrix
      datamatrix[[peak]]<-newSampleMatrix
    }
    lobject<-methods::new("LDO", name=name, dataMatrices=datamatrix, backgroundMatrices=bgmatrix, peaknames=peaknames, times=times, ids=ids, labels=as.factor(labels));
  }else{
    datamatrix<-list();
    peak<-peaknames[1]
    for(peak in peaknames){
      newSampleMatrix <- data.frame(class=factor(longData[,class], levels=labels), 
                                    id = factor(longData[,id], levels=ids), 
                                    time=as.numeric(longData[,time]), 
                                    value=as.numeric(longData[,peak]))
      datamatrix[[peak]]<-newSampleMatrix
    }
    lobject<-methods::new("LDO", name=name, dataMatrices=datamatrix, peaknames=peaknames, times=times, ids=ids, labels=as.factor(labels));
  }
  return(lobject);
}


#' Create a new 'LDO' Object that only contains the selected components.
#' 
#' @param ldo Longitudinal Data Object
#' @param components Component names to select for the new ldo object. Only elements from this list that overlap with the peak names in the given ldo, are utilized.
#' @param name Name of newly created 'LDO' object.
#' @importFrom methods new
#' @export
#' @return new ldo object only containing the selected components.
#' 
selectComponents=function(ldo, components, name=paste(ldo@name, ' selected')){
  
  if(length(components)<1){
    stop("The components list has to contain at least one name");
  }
  
  peaknames= ldo@peaknames[ldo@peaknames %in% (components)]
  
  if(sum(!(components %in% peaknames))>0){
    warning("Some of the components are not available in this ldo!");
  }
  
  # Select the components and adjust the data and background matrix.
  datamatrix = ldo@dataMatrices[components]
  bgmatrix = ldo@backgroundMatrices[components]
  
  # Create new LDO Object with adjusted components.
  lobject<-methods::new("LDO", name=name, dataMatrices=datamatrix,backgroundMatrices=bgmatrix, peaknames=components, times=ldo@times, ids=ldo@ids, labels=ldo@labels);

  return(lobject);
}


#' Get colors for the plotting function.
#' 
#' @param label class labels of the samples
#' @param size size of the color vector to be created
#' @return col vector of colors created
#' 
getColor=function(label, size){
  # size<-10
  if(label==1){
    col<- grDevices::colorRampPalette(c( "darkgoldenrod3", "yellow", "khaki"))(size);
  }else   if(label==2){
    col<- grDevices::colorRampPalette(c("darkblue", "cadetblue1"))(size);
  }else   if(label==3){
    col<- grDevices::colorRampPalette(c("darkgreen", "green", "lightgreen"))(size);
  }else {
    col<- grDevices::colorRampPalette(c("darkred", "red", "pink"))(size);
  }
  return(col)
}