R/tidy_exwas_outcomes.R

In groovearmada/exwas:

tidy_exwas_outcomes <- function(input_data, variable_list, exwas_result, subcohort, file_name) {

  # # Transform the results so they can be used for plotting and
  # # exported as .csv table

  #calculate effective number of
  # tests (Meff), to calculate FDR corrected p value)
  if (class(input_data) == "mids") {
    cormat <- mice::complete(input_data) %>%
      dplyr::select(variable_list$Term) %>%

      # calculate correlation matrix
      cor(use = "pairwise.complete.obs")

  } else {

    # calculate correlation matrix
    cormat <- dplyr::select(input_data, ends_with("_log2_iqr"), ends_with("_ln_iqr")) %>%
      cor(use = "pairwise.complete.obs")
  }

  # meff calculation based on Li 2012 method
  lambdas <- base::eigen(cormat)$values
  Meff <- ncol(cormat) - sum((lambdas > 1) * (lambdas - 1))  #32.69156

  # Add p value corrected for multiple comparison for
  # correlated explanatory variables, Li 2012 method
  exwas_result <- exwas_result %>%
    dplyr::mutate(p_value_FDR = pmin(1, Meff * p_value),
                  p_value_round = round(p_value, 3))

  if (is.null(subcohort) == TRUE) {

    # Add exposure names
    tidy_exwas_res <- exwas_result %>%
      dplyr::mutate(Exposure = variable_list[, 2]) %>%

      # add variable descriptions
      dplyr::left_join(variable_list[, c(2:3)], ., by = "Exposure") %>%

      # Order the table by Exposure
      dplyr::arrange(Family, Exposure) %>%

      # change the order of variables so it is ready for the table
      .[c(2:47, 53, 48:52, 1), ]

    # remove unnecessary variables
    pooled_fits_plot <- dplyr::select(tidy_exwas_res, Exposure,
                                      Estimate, conf_low, conf_high, p_value)

    exwas_sel <- dplyr::select(tidy_exwas_res, Family, Exposure, estCI,
                               p_value_round, p_value_FDR)

  } else {

    # add variable descriptions
    lev_num <- dplyr::select(exwas_result, sub_cohort) %>%
      unique() %>%
      nrow()

    tidy_exwas_res <- exwas_result %>%
      dplyr::mutate(Exposure = rep(variable_list[, 2], times = lev_num)) %>%
      dplyr::left_join(variable_list[, c(2:3)], ., by = "Exposure") %>%

      # Order the table by subcohort and Exposure
      dplyr::arrange(sub_cohort, Family, Exposure)

    pooled_fits_plot <- dplyr::select(tidy_exwas_res, sub_cohort,
                                      Exposure, Estimate, conf_low, conf_high, p_value)

    # Order the table by subcohort and Exposure
    exwas_sel <- dplyr::select(tidy_exwas_res, sub_cohort, Family, Exposure,
                               estCI, p_value_round, p_value_FDR)
  }

  # export table result to a .csv file
  write.csv(pooled_fits_plot, file = paste0(file_name, "_plot.csv"))

  # export table result to a .csv file
  write.csv(exwas_sel, file = paste0(file_name, ".csv"))

  # return table ready to plot
  return(pooled_fits_plot)
}