R/StepThree.R

In SMAHP: Survival Mediation Analysis of High-Dimensional Proteogenomic Data

#' Post-hoc model
#'
#' An auxiliary function which generates p-value based on outcome and mediation model.
#' @param X An n by p matrix of exposures.
#' @param M An n by p matrix of mediators.
#' @param C An n by p matrix of covariates. If there are no covariates, set C = NULL.
#' @param time A vector of survival time of samples.
#' @param status A vector of status indicator: 0=alive, 1=dead.
#' @param X_sel_Y_s1 Outputs from StepOne: A vector of selected X for Y.
#' @param M_X_sel_s2 Outputs from StepTwo: A data table with selected M, X pairs and related effect size.
#' @return A list with the following components:
#'   \item{p_beta_m}{p-values generated from outcome model}
#'   \item{p_alpha_x}{p-values generated from mediation modell}
#'   \item{outcome_model}{coefficient estimates from outcome model}
#'   \item{med_results}{coefficient estimates from mediation model}
#' @import survival
#' @export
#' @examples
#' \donttest{
#' data(example_dat)
#' surv_dat <- example_dat$surv_dat
#' res_step1 <- StepOne(X = example_dat$X, M = example_dat$M, time = surv_dat$time,
#' status = surv_dat$status, model_option = "MCP")
#' M_X_sel_s2 <- StepTwo(X = example_dat$X, M = example_dat$M, time = surv_dat$time,
#' status = surv_dat$status, X_sel_Y_s1 = res_step1$X_sel_Y_s1, M_X_s1 = res_step1$M_X_s1,
#' M_sel_Y_s1 = res_step1$M_sel_Y_s1)
#' res_step3 <- StepThree(X = example_dat$X, M = example_dat$M, C = example_dat$C,
#' time = surv_dat$time, status = surv_dat$status, X_sel_Y_s1 = res_step1$X_sel_Y_s1,
#' M_X_sel_s2 = M_X_sel_s2)
#' }

StepThree <- function(X, M, C, time, status, X_sel_Y_s1, M_X_sel_s2){
  X_u <- X
  M_u <- M
  colnames(X_u) <- paste0("X_", colnames(X))
  colnames(M_u) <- paste0("M_", colnames(M))

  d_x <- ncol(X)
  M_sel_u <- paste0("M_", unique(M_X_sel_s2$M))
  X_sel_Y_s1 <- paste0("X_", X_sel_Y_s1)
  n_m_s3 <- length(M_sel_u)
  surv_data <- data.frame(surv_time = time, status = status)
  surv_model_data_s3 <- cbind(M_u[ ,M_sel_u], X_u[ ,X_sel_Y_s1], C, surv_data)
  surv_model_s3 <- survreg(Surv(surv_time, status) ~ ., data = surv_model_data_s3, dist = "lognormal")
  mod_results <- surv_model_s3$coefficients
  # transform p
  M_pos <- which(startsWith(names(mod_results), "M"))
  if(length(M_sel_u) == 1){
    M_pos <- 2
  }
  p_list <- abs(summary(surv_model_s3)$table[M_pos,1])/summary(surv_model_s3)$table[M_pos,2]
  p_list <- 2*(1 - pnorm(p_list))
  p_beta_m <- p_list

  col_X <- colnames(X)
  M_sel <- unique(M_X_sel_s2$M)
  p_matrix <- matrix(nrow = d_x, ncol = length(M_sel))
  med_results <- list()
  for(i in 1:length(M_sel)){
    M_k <- M_sel[i]
    X_sel_M_k_s3 <- M_X_sel_s2[M == M_k, X]
    lm_model_data_s3 <- cbind(M[ ,M_k], X[ ,X_sel_M_k_s3], C)
    lm_model_data_s3 <- as.data.table(lm_model_data_s3)
    lm_model_s3 <- lm(V1 ~ ., data = lm_model_data_s3)
    lm_model_s3_results <- summary(lm_model_s3)$coefficients[ ,1]

    med_results[[M_k]] <- lm_model_s3_results

    if(is.null(C)){
      n_c = 0
    } else {
      n_c <- ncol(C)
    }
    X_pos <- c(3:(length(lm_model_s3_results) - n_c))
    X_name <- names(lm_model_s3_results)[X_pos]

    p_list <- abs(summary(lm_model_s3)$coefficients[X_pos,1])/summary(lm_model_s3)$coefficients[X_pos, 2]
    p_list <- 2*(1 - pnorm(p_list))
    p_list <- as.matrix(p_list)
    if(length(X_name == 1)){
      rownames(p_list) <- X_name
    }
    col_p_list_r <- setdiff(col_X, X_name)
    p_list_r <- matrix(NA, nrow = length(col_p_list_r), ncol = 1)
    rownames(p_list_r) <- col_p_list_r
    p_list_final <- rbind(p_list, p_list_r)
    p_list_final <- p_list_final[colnames(X), ]
    p_matrix[ ,i] <- p_list_final
  }
  rownames(p_matrix) <- colnames(X)
  colnames(p_matrix) <- M_sel
  p_alpha_x <- p_matrix

  return(list(p_beta_m = p_beta_m, p_alpha_x = p_alpha_x, outcome_model = mod_results, med_results = med_results))
}