add wald-test for glm

src/v2/HTRSIM/R/modelFitting.R

+#' Launch wald test
+#'
+#' @param model.res output of glm.nb
+#' @param term index of the coeff to test
+#' @param threshold HO : coeff < threshold
+#' @param initvec vector of 0
+#' @import glmglrt
+#' @return a numeric pvalue
+#' @export
+#'
+#' @examples
+test_wald <- function(model.res, term, threshold, initvec){
+  constrast_vec = replace(initvec, term , 1)
+  model.res$coefficients <- abs(model.res$coefficients)
+  wald.test.pvalue = glmglrt::p_value_contrast(model.res, contrast = constrast_vec ,  alternative = "greater", method = "Wald", H0 = threshold)
+  return(wald.test.pvalue %>% as.numeric())
+}
+
+
+
 #' Launch MASS:GLM.NB
 #'
 #' @param gene_count a row of kij simulated table
@@ -5,22 +25,36 @@
 #' @import dplyr
 #' @import stringr
 #' @import stats
+#' @import purrr
+#' @import furrr
 #' @return a dtf
 #' @export
 #'
 #' @examples
-reshapeGlmRes <- function(fit, i, error_bool = F){
+reshapeGlmRes <- function(fit, i, threshold  ,  error_bool = F){
     if (error_bool == T){ ## error while fiting model
       res = list(Inference = NA, pval = NA,
-                 beta = NA, gene_id = paste("gene", i, sep = ""), type = NA, deviance = NA) %>%
+                 beta = NA, gene_id = paste("gene", i, sep = ""), type = NA, deviance = NA, p.val = NA) %>%
             data.frame()
     }
     else{ ## success to fit model
+
+      ### wald test ###
+      message(paste("WALD test:\n HO:", '|Beta| <' , threshold, sep = ' ' ))
+      vecofzero = rep(0, length(coef(fit)))
+      list_pval = 1:length(coef(fit)) %>% purrr::map(.x = .,  ~test_wald(model.res = fit, term= .x ,  threshold , initvec = vecofzero) ) %>% unlist()
+
+
+
+      ###
       res = coef(summary(fit))[,c(1,4)] %>%
               data.frame() %>%
-              dplyr::rename(., pval = "Pr...z.." , Inference = "Estimate") %>%
+              dplyr::rename(.,  Inference = "Estimate") %>%
               dplyr::mutate(beta = stringr::str_remove_all(rownames(.), "[//(//)]")) %>%
-              dplyr::mutate(beta = stringr::str_replace(beta, ":", "."))
+              dplyr::mutate(beta = stringr::str_replace(beta, ":", ".")) %>%
+              dplyr::mutate(pval = list_pval) %>%
+              dplyr::select(-"Pr...z..")
+
       rownames(res) <- NULL
       res = res  %>%
               dplyr::mutate(gene_id = paste("gene", i, sep = "") ) %>%
@@ -48,21 +82,24 @@ reshapeGlmRes <- function(fit, i, error_bool = F){
 #' @export
 #'
 #' @examples
-run.glm <- function(gene_count, i, design.dtf) {
+run.glm <- function(gene_count, i, design.dtf, threshold = 0) {
   y = gene_count
   genotype = design2simulate$design2simulate$genotype
   environment = design2simulate$design2simulate$environment
+
+  message("Fitting model ...")
+
   df_gene_i = list(y = y , genotype = genotype,environment = environment) %>% data.frame()
   rownames(df_gene_i) <- NULL
-  print(i)
+  #print(i)
   tryCatch({
 
     fit = MASS::glm.nb(y ~ genotype + environment + genotype:environment, data = df_gene_i, link = log)
-    return(reshapeGlmRes(fit, i))
+    return(reshapeGlmRes(fit, i, threshold))
   },
 
   error = function(cnd){
-    return(reshapeGlmRes(NULL, i, error_bool = T))
+    return(reshapeGlmRes(NULL, i, threshold  , error_bool = T))
 }
 )
 }
@@ -86,3 +123,7 @@ run.deseq <- function(tabl_cnts, bioDesign, model = ~ genotype + environment + g
   return(dds)
 
 }
+
+
+
+