init HTRfit

src/model.fit/HTRfit/R/model_fitting.R

+
+#' launch glm on data
+#' @param data2fit
+#' @param model2fit
+#' @param fit_by
+#' @param threads
+#' @import future
+#' @import furrr
+#' @return fit list 
+#' @export
+#'
+#' @examples
+launch.glm <- function( data2fit, model2fit = k_ij ~ genotype + environment + genotype:environment, 
+                    fit_by = "gene_id", threads = 4 ) {
+
+    iteration_list = data2fit[fit_by] %>% unique() %>% as.vector() %>% unname() %>% .[[1]]  
+    future::plan(multisession, workers = threads)
+    fit_list = iteration_list %>% furrr::future_map(.x = ., 
+                                ~fit.glm( data2fit[data2fit[fit_by] == .x,], 
+                                          .x, 
+                                           model2fit) )
+    return(fit_list)
+
+}
+
+
+#' fit model on data
+#' @param data 
+#' @param id
+#' @import MASS
+#' @import dplyr
+#' @return fit 
+#' @export
+#'
+#' @examples
+fit.glm <- function(data, id, model2fit) {
+  fit = MASS::glm.nb(model2fit, data = data, link = log)
+  fit.dtf = tidySummary(fit)
+  fit.dtf$inference = fit.dtf$inference %>% dplyr::mutate(gene_id = id)
+  fit.dtf$fitQuality = fit.dtf$fitQuality %>% dplyr::mutate(gene_id = id)
+  return(fit.dtf)
+}
+
+#' fit model on data
+#' @param fit 
+#' @import broom
+#' @return list of element 
+#' @export
+#'
+#' @examples
+tidySummary <- function(fit){
+    dtf = broom::tidy(fit)
+    gl = broom::glance(fit)
+    return( list(inference = dtf, fitQuality = gl) )
+}
+
+#' convert list to dataframe
+#' @param fit 
+#' @return list of dtf
+#' @export
+#'
+#' @examples
+listFit2dtf <- function(list_fit){
+    tmp = do.call(cbind, list_fit)
+    inference.dtf = do.call(rbind, tmp[1,])
+    fitQuality.dtf = do.call(rbind, tmp[2,])
+    return( list(inference = inference.dtf, fitQuality = fitQuality.dtf) )
+}
\ No newline at end of file

src/model.fit/HTRfit/R/preprocessing.R

+#' get df2fit for a given gene
+#' @param k_ij 
+#' @param design_simulation
+#' @param gene_name
+#' @import dplyr
+#' @return dataframe 
+#' @export
+#'
+#' @examples
+getdf2fit <- function(k_ij , design_simulation, gene_name){
+  df_gene_i = cbind(design_simulation, k_ij) %>% 
+                        dplyr::mutate(gene_id = gene_name) %>% 
+                        dplyr::mutate(gene_id = gene_name)
+  rownames(df_gene_i) <- NULL
+  return(df_gene_i)
+}   
+
+#' get df2fit for a given gene
+#' @param countTable 
+#' @param experimental_design
+#' @import purrr
+#' @return dataframe 
+#' @export
+#'
+#' @examples
+reshapeCounTable<- function(countTable, experimental_design ){
+
+    gene_id_list = rownames( countTable )
+    list_df2fit = purrr::map( .x = gene_id_list, ~getdf2fit( countTable[ .x, ], 
+                                           experimental_design, .x) )
+    data2fit = do.call( rbind, list_df2fit )
+    return(data2fit)
+}
\ No newline at end of file

src/model.fit/HTRfit/R/wald_test.R

+#' fit model on data
+#' @param W 
+#' @param threshold
+#' @param alternative
+#' @import stats
+#' @return pvalue 
+#' @export
+#'
+#' @examples
+wald_test <- function( W, threshold = 0, alternative = "two.sided"){
+    ## two.sided : P ( X > threshold ) + P ( X < threshold )
+    pval = 2*(stats::pnorm( W, mean = 0, sd = 1, lower.tail = FALSE))
+
+    ## greater : P ( X > threshold )
+    pval = stats::pnorm( W, mean = 0, sd = 1, lower.tail = FALSE)
+    
+    ## less : P ( X <= threshold )
+    pval = stats::pnorm( W , mean = 0, sd = 1, lower.tail = TRUE)
+    return(pval)
+}
+
+
+#' fit model on data
+#' @param estimate 
+#' @param threshold
+#' @param stdError
+#' @return W 
+#' @export
+#'
+#' @examples
+getStatisticWaldTest<- function(estimate,stdError, threshold ){
+    ## cf https://en.wikipedia.org/wiki/Wald_test
+    return( ( estimate - threshold ) / stdError )
+}
+
+
+

src/model.fit/HTRfit/devtools_history.R

+usethis::use_build_ignore("devtools_history.R")
+usethis::use_package('tidyverse', type = "depends")
+usethis::use_package("dplyr")
+usethis::use_package("purrr")
+usethis::use_package("furrr")
+usethis::use_package("MASS")
+usethis::use_package("future")
+usethis::use_package("stats")
+usethis::use_package("broom")
+
+
+
+