prolfqua: Proteomics Label Free Quantification

# Contrasts -----

#' Estimate contrasts using Wald Test
#' @export
#' @family modelling
#' @examples
#'
#' # Fitting mixed effects model to peptide data
#' istar <- prolfqua::sim_lfq_data_peptide_config()
#'
#' modelFunction <-
#' strategy_lmer("abundance  ~ group_ + (1 | peptide_Id) + (1 | sample)")
#'
#' config <- istar$config
#' config$table$hierarchy_keys_depth()
#'
#' mod <- build_model(
#'  istar$data,
#'  modelFunction,
#'  subject_Id = config$table$hierarchy_keys_depth()
#'  )
#'
#' prolfqua::model_summary(mod)
#'  Contr <- c("groupA_vs_Ctrl" = "group_A - group_Ctrl",
#'     "dil.e_vs_b" = "group_A - group_Ctrl" )
#' #prolfqua::Contrasts$debug("get_contrasts")
#' contrastX <- prolfqua::Contrasts$new(mod, Contr)
#' contrastX$get_linfct(global=FALSE)
#' contrastX$get_contrasts()
#' contrastX$get_contrast_sides()
#' contrastX$column_description()
#'
Contrasts <- R6::R6Class(
  "Contrast",
  inherit = ContrastsInterface,
  public = list(
    #' @field models Model
    models = NULL,
    #' @field contrasts character with contrasts
    contrasts = character(),
    #' @field contrastfun function to compute contrasts
    contrastfun = NULL,
    #' @field modelName model name
    modelName = character(),
    #' @field subject_Id name of column containing e.g., protein Id's
    subject_Id = character(),
    #' @field p.adjust function to adjust p-values (default prolfqua::adjust_p_values)
    p.adjust = NULL,
    #' @field contrast_result data frame containing results of contrast computation
    contrast_result = NULL,
    #' @field global use a global linear function (determined by get_linfct)
    global = TRUE,
    #' @field protein_annot holds protein annotation
    protein_annot = NULL,
    #' @description
    #' initialize
    #' create Contrast
    #' @param model a dataframe with a structure similar to that generated by \code{\link{build_model}}
    #' @param contrasts a character vector with contrast specificiation
    #' @param p.adjust function to adjust the p-values
    #' @param global development/internal argument (if FALSE determine linfct for each model.)
    #' @param modelName name of contrast method, default WaldTest
    initialize = function(model,
                          contrasts,
                          p.adjust = prolfqua::adjust_p_values,
                          global = FALSE,
                          modelName = "WaldTest"
    ){
      self$models = model$modelDF |> dplyr::filter(exists_lmer == TRUE)
      self$contrasts = contrasts
      self$contrastfun = model$model_strategy$contrast_fun
      self$modelName =  modelName
      self$subject_Id = model$subject_Id
      self$p.adjust = p.adjust
      self$global = global
    },
    #' @description
    #' get both sides of contrasts
    get_contrast_sides = function(){
      # extract contrast sides
      tt <- self$contrasts[grep("-",self$contrasts)]
      tt <- tibble(contrast = names(tt) , rhs = tt)
      tt <- tt |> mutate(rhs = gsub("[` ]","",rhs)) |>
        tidyr::separate(rhs, c("group_1", "group_2"), sep = "-")
      return(tt)
    },
    #' @description
    #' get linear functions from contrasts
    #' @param global logical TRUE - get the a linear functions for all models, FALSE - linear function for each model
    get_linfct = function(global = TRUE){
      linfct <- function(model, contrast){
        linfct <- linfct_from_model(model, as_list = FALSE)
        linfct <- unique(linfct) # needed for single factor models
        # namtmp <- paste0("avg_",names(self$contrasts))
        # tmp <- paste0("(", gsub(" - ", " + ", self$contrasts), ")/2")
        # names(tmp) <- namtmp
        # cntrasts <- c(self$contrasts, tmp)
        linfct_A <- linfct_matrix_contrasts(linfct, contrast)
        return( linfct_A )
      }
      if (global) {

        model <- get_complete_model_fit(self$models)$linear_model[[1]]
        res <- linfct( model, self$contrasts )
        return( res )

      }else{

        res <- vector(mode = "list", nrow(self$models))
        for (i in seq_along(self$models$linear_model)) {
          res[[i]] <- linfct(self$models$linear_model[[i]], contrast = self$contrasts)
        }
        return(res)
      }
    },
    #' @description
    #' get table with contrast estimates
    #' @param all should all columns be returned (default FALSE)
    #' @return data.frame with contrasts
    #'
    get_contrasts = function(all = FALSE){

      if (is.null(self$contrast_result) ) {
        message("determine linear functions:")
        linfct <- self$get_linfct(global = self$global)
        #contrast_sides <- self$get_contrast_sides()
        message("compute contrasts:")
        # TODO (goes into calling code)
        contrast_result <- contrasts_linfct(
          self$models,
          linfct,
          subject_Id = self$subject_Id,
          contrastfun = self$contrastfun )
        contrast_result <- ungroup(contrast_result)

        contrast_result <- dplyr::rename(contrast_result, contrast = lhs, diff = estimate)

        differences <- contrast_result |>
          dplyr::filter(contrast %in% names(self$contrasts))

        avgAbd <- contrast_result |>
          dplyr::select(dplyr::all_of(c(self$subject_Id, "contrast", "diff"))) |>
          dplyr::filter(startsWith(contrast , "avg_"))

        avgAbd$contrast <- gsub("^avg_","", avgAbd$contrast)
        avgAbd <- avgAbd |> dplyr::rename(avgAbd = diff)
        contrast_result <- left_join(differences, avgAbd)

        contrast_result <- self$p.adjust(contrast_result,
                                         column = "p.value",
                                         group_by_col = "contrast")
        contrast_result <- contrast_result |> relocate("FDR", .after="diff")
        contrast_result <- mutate(contrast_result, modelName = self$modelName, .before = 1)
        self$contrast_result <- contrast_result

      }
      res <- if (!all) {
        self$contrast_result |>
          select( -all_of(c("sigma.model",
                            "df.residual.model",
                            "isSingular")))
      }else{
        self$contrast_result
      }

      stopifnot(all(super$column_description()$column_name %in% colnames(res)))
      return(res)

    },
    #' @description
    #' return \code{\link{ContrastsPlotter}}
    #' creates Contrast_Plotter
    #' @param FCthreshold fold change threshold to show in plots
    #' @param FDRthreshold FDR threshold to show in plots
    #' @return \code{\link{ContrastsPlotter}}
    get_Plotter = function(
    FCthreshold = 1,
    FDRthreshold = 0.1
    ){
      contrast_result <- self$get_contrasts()
      res <- ContrastsPlotter$new(
        contrast_result,
        subject_Id = self$subject_Id,
        fcthresh = FCthreshold,
        volcano = list(list(score = "p.value", thresh = FDRthreshold),
                       list(score = "FDR", thresh = FDRthreshold)),
        histogram = list(list(score = "p.value", xlim = c(0,1,0.05)),
                         list(score = "FDR", xlim = c(0,1,0.05))),
        score = list(list(score = "statistic", thresh = 5)),
        modelName = "modelName",
        diff = "diff",
        contrast = "contrast"
      )
      return(res)
    },
    #' @description
    #' convert to wide format
    #' @param columns value column default p.value
    #' @return data.frame
    to_wide = function(columns = c("p.value", "FDR", "statistic")){
      contrast_minimal <- self$get_contrasts()
      contrasts_wide <- pivot_model_contrasts_2_Wide(contrast_minimal,
                                                     subject_Id = self$subject_Id,
                                                     columns = c("diff", columns),
                                                     contrast = 'contrast')
      return(contrasts_wide)
    }
  ))