prolfqua: Proteomics Label Free Quantification

#'
#'
#'
.linfct <- function(model, contrast, avg = TRUE){
  linfct <- linfct_from_model(model, as_list = FALSE)
  linfct <- unique(linfct) # needed for single factor models
  if (avg) {
    namtmp <- paste0("avg_",names(contrast))
    cntr_avg <- paste0("(", gsub(" - ", " + ", contrast), ")/2")
    names(cntr_avg) <- namtmp
    contrast <- c(contrast, cntr_avg)
  }
  linfct_A <- linfct_matrix_contrasts(linfct, contrast)
  return( linfct_A )
}

# Contrasts -----

#' Estimate contrasts using Wald Test
#' @export
#' @family modelling
#' @examples
#'
#' modi <- sim_build_models_logistf(model = "parallel3", weight_missing = 1)
#' contrasts <- c(Avs = "group_A - group_B", AvsCtrl = "group_A - group_Ctrl")
#'
#' ctr <- ContrastsFirth$new(modi,contrasts)
#' ctr$get_contrast_sides()
#' ctr$get_linfct()
#' ctr$get_contrasts()
#'
#' mod3 <- sim_build_models_logistf(model = "parallel3", weight_missing = 1, peptide=TRUE)
#'
#' # ContrastsFirth$debug("get_contrasts")
#' ctrpep <- ContrastsFirth$new(mod3,contrasts)
#' ctrpep$get_contrast_sides()
#'
#' ctrpep$get_linfct()
#' ctrpep$get_contrasts()
#' pl <- ctrpep$get_Plotter()
#'
ContrastsFirth <- R6::R6Class(
  "ContrastFrith",
  inherit = ContrastsInterface,
  public = list(
    #' @field models Model
    models = NULL,
    #' @field contrasts character with contrasts
    contrasts = character(),
    #' @field contrastfun function to compute contrasts
    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,
    #' @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 modelName name of contrast method, default WaldTest
    initialize = function(model,
                          contrasts,
                          p.adjust = prolfqua::adjust_p_values,
                          modelName = "WaldTestFirth"
    ){
      self$models = model
      self$contrasts = contrasts
      self$modelName =  modelName
      self$subject_Id = model$subject_Id
      self$p.adjust = p.adjust
    },
    #' @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 avg logical TRUE - get also linfct for averages
    get_linfct = function(avg = TRUE){

      if (!is.null(self$models$models$models1)) {
        model1DF <- self$models$models$models1$modelDF
        res1 <- vector(mode = "list", nrow(model1DF))
        pb <- progress::progress_bar$new(total = length(model1DF$linear_model))
        model <- get_complete_model_fit(model1DF)$linear_model[[1]]
        compmodel <- .linfct( model, self$contrasts, avg = avg )
        mcoef <- max(model1DF$nrcoeff_not_NA)

        for (i in seq_along(model1DF$linear_model)) {
          pb$tick()
          res1[[i]] <-  if (model1DF$nrcoeff_not_NA[[i]] == mcoef) { compmodel } else {
            .linfct(model1DF$linear_model[[i]],
                    contrast = self$contrasts, avg = avg)}
        }
        self$models$models$models1$modelDF$linfct <- res1
      }

      if (!is.null(self$models$models$models2)) {
        model2DF <- self$models$models$models2$modelDF
        pb <- progress::progress_bar$new(total = length(model2DF$linear_model))
        res2 <- vector(mode = "list", nrow(model2DF))
        for (i in seq_along(model2DF$linear_model)) {
          pb$tick()
          res2[[i]] <-
            .linfct(model2DF$linear_model[[i]],
                    contrast = self$contrasts, avg = avg)
        }
        self$models$models$models2$modelDF$linfct <- res2
      }
      return(self)
    },
    #' @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:")
        self$get_linfct()
        message("get_contrasts -> contrasts_linfct")
        # TODO (goes into calling code)
        models1 <- self$models$models$models1
        contrast_result1 <- NULL
        contrast_result2 <- NULL
        if (!is.null(models1)) {
          contrast_result1 <- contrasts_linfct_firth(models1)
          contrast_result1 <- ungroup(contrast_result1)
        }
        models2 <- self$models$models$models2
        if (!is.null(models2)) {
          contrast_result2 <- contrasts_linfct_firth(models2)
          contrast_result2 <- ungroup(contrast_result2)
        }
        contrast_result <- bind_rows(contrast_result2, contrast_result1)
        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)
    }
  ))