diff --git a/NAMESPACE b/NAMESPACE
index bd1c8e33dc27762287f9d98187f4add7f2503186..6b45d02d6500eb66090777e854c3640c520bd303 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -95,6 +95,7 @@ export(is_formula_mixedTypeI)
 export(is_fullrank)
 export(is_mixedEffect_inFormula)
 export(is_positive_definite)
+export(is_validGroupBy)
 export(join_dtf)
 export(launchFit)
 export(launchUpdate)
@@ -105,6 +106,7 @@ export(parallel_fit)
 export(parallel_update)
 export(prepareData2computeInteraction)
 export(prepareData2fit)
+export(prepare_dataParallel)
 export(removeDigitsAtEnd)
 export(removeDuplicatedWord)
 export(renameColumns)
@@ -118,7 +120,6 @@ export(scaleCountsTable)
 export(set_correlation)
 export(simulationReport)
 export(subsetByTermLabel)
-export(subsetData_andfit)
 export(subsetFixEffectInferred)
 export(subsetGenes)
 export(subset_glance)
diff --git a/R/fitmodel.R b/R/fitmodel.R
index b38b1df102da2b147a3ba260d89d2c791439a288..74e1f5105beed8abfb771fee91582557bb81c01c 100644
--- a/R/fitmodel.R
+++ b/R/fitmodel.R
@@ -25,6 +25,27 @@ isValidInput2fit <- function(data2fit, formula){
 }
 
 
+
+#' Check if group by exist in data
+#'
+#' @param data The data framecontaining the variables to be used for model fitting.
+#' @param group_by Column name in data representing the grouping variable 
+#'
+#' @return \code{TRUE} if exist otherwise an error is raised
+#'
+#' @examples
+#' is_validGroupBy(mtcars, 'mpg')
+#' @export
+is_validGroupBy <- function(data, group_by){
+  validGroupBy <- group_by %in% names(data)
+  if (!validGroupBy) 
+    stop("<Group by> doen't exist in data !")
+  return(TRUE)
+}
+
+
+
+
 #' Drop Random Effects from a Formula
 #'
 #' This function allows you to remove random effects from a formula by specifying 
@@ -105,7 +126,7 @@ is_fullrank <- function(metadata, formula) {
 
 
 #' Fit a model using the fitModel function.
-#' @param group group id to save in glmmTMB obj (usefull for update !)
+#' @param group ID to fit
 #' @param formula Formula specifying the model formula
 #' @param data Data frame containing the data
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
@@ -113,7 +134,7 @@ is_fullrank <- function(metadata, formula) {
 #' @export
 #' @examples
 #' fitModel("mtcars" , formula = mpg ~ cyl + disp, data = mtcars)
-fitModel <- function(group, formula, data, ...) {
+fitModel <- function(group , formula, data, ...) {
   # Fit the model using glm.nb from the GLmmTMB package
   model <- glmmTMB::glmmTMB(formula, ..., data = data ) 
   model$frame <- data
@@ -132,24 +153,22 @@ fitModel <- function(group, formula, data, ...) {
 
 #' Fit the model based using fitModel functions.
 #'
-#' @param group The specific group to fit the model for
+#' @param groups list of group ID
 #' @param group_by Column name in data representing the grouping variable
-#' @param formula Formula specifying the model formula
 #' @param data Data frame containing the data
-#' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
-#' @return Fitted model object or NULL if there was an error
+#' @return list of dataframe 
 #' @export
+#' @importFrom stats setNames
 #' @examples
-#' subsetData_andfit(group = "setosa", group_by = "Species", 
-#'                  formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
+#' prepare_dataParallel(groups = iris$Species, group_by = "Species", 
 #'                  data = iris )
-subsetData_andfit <- function(group, group_by, formula, data, ...) {
-  subset_data <- data[data[[group_by]] == group, ]
-  fit_res <- fitModel(group, formula, subset_data, ...)
-  #glance_df <- glance.negbin(group_by ,group , fit_res)
-  #tidy_df <- tidy.negbin(group_by ,group,fit_res )
-  #list(glance = glance_df, summary = tidy_df)
-  fit_res
+prepare_dataParallel <- function(groups, group_by, data) {
+  
+  l_data2parallel <- lapply( stats::setNames(groups, groups) , function( group_id ){
+                      subset_data <- data[ data[[ group_by ]] == group_id, ]
+                      return(subset_data)
+                  })
+  return(l_data2parallel)
 }
 
 
@@ -159,22 +178,22 @@ subsetData_andfit <- function(group, group_by, formula, data, ...) {
 #' This function fits the model using the specified group, group_by, formula, and data.
 #' It handles warnings and errors during the fitting process and returns the fitted model or NULL if there was an error.
 #'
-#' @param group The specific group to fit the model for
+#' @param data Data frame containing the data
 #' @param group_by Column name in data representing the grouping variable
 #' @param formula Formula specifying the model formula
-#' @param data Data frame containing the data
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
 #' @return List with 'glance' and 'summary' attributes representing the fitted model or NULL if there was an error
 #' @export
 #' @examples
-#' launchFit(group = "setosa", group_by = "Species", 
+#' launchFit(group_by = "Species", 
 #'            formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
-#'            data = iris )
-launchFit <- function(group, group_by, formula, data, ...) {
+#'            data = iris[ iris[["Species"]] == "setosa" , ] )
+launchFit <- function(data, group_by, formula, ...) {
+  group <- unique(data[[ group_by ]]) 
   tryCatch(
     expr = {
       withCallingHandlers(
-          subsetData_andfit(group, group_by, formula, data, ...),
+          fitModel(group , formula, data, ...),
           warning = function(w) {
             message(paste(Sys.time(), "warning for group", group, ":", conditionMessage(w)))
             invokeRestart("muffleWarning")
@@ -183,7 +202,6 @@ launchFit <- function(group, group_by, formula, data, ...) {
     error = function(e) {
       message(paste(Sys.time(), "error for group", group, ":", conditionMessage(e)))
       NULL
-      #return(list(glance = empty.glance.negbin(group_by, group), summary = empty.tidy.negbin(group_by, group)))
     }
   )
 }
@@ -201,10 +219,9 @@ launchFit <- function(group, group_by, formula, data, ...) {
 #' @param log_file File to write log (default : Rtmpdir/htrfit.log)
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
 #' @return List of fitted model objects or NULL for any errors
-#' @importFrom stats setNames
 #' @export
 #' @examples
-#' parallel_fit(group_by = "Species", "setosa", 
+#' parallel_fit(group_by = "Species", groups =  iris$Species, 
 #'                formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
 #'                data = iris, n.cores = 1 )
 parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL, 
@@ -212,10 +229,14 @@ parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL,
   
   if (is.null(n.cores)) n.cores <- parallel::detectCores()
   
+  ## get data for parallelization
+  l_data2parallel <- prepare_dataParallel(groups, group_by, data)
+
   clust <- parallel::makeCluster(n.cores, outfile = log_file)
-  parallel::clusterExport(clust, c("subsetData_andfit", "fitModel"),  envir=environment())
-  results_fit <- parallel::parLapply(clust, X = stats::setNames(groups, groups), fun = launchFit, 
-                      group_by = group_by, formula = formula, data = data, ...)
+  parallel::clusterExport(clust, c("fitModel"),  envir=environment())
+  results_fit <- parallel::parLapply(clust, X = l_data2parallel, 
+                                     fun = launchFit, 
+                                     group_by = group_by, formula = formula, ...)
                                      
   parallel::stopCluster(clust)
   #closeAllConnections()
@@ -239,15 +260,16 @@ parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL,
 #'                  data = iris, group_by = "Species", n.cores = 1) 
 fitModelParallel <- function(formula, data, group_by, n.cores = NULL, log_file = paste(tempdir(check = FALSE), "htrfit.log", sep = "/"), ...) {
   
-  ## SOme verification
+  ## Some verification
   isValidInput2fit(data, formula)
   is_fullrank(data, formula)
-  
+  is_validGroupBy(data, group_by)
+
   ## -- print log location
   message( paste("Log file location", log_file, sep =': ') ) 
   
-  groups <- unique(data[[group_by]])
   # Fit models in parallel and capture the results
+  groups <- unique(data[[ group_by ]])
   results <- parallel_fit(groups, group_by, formula, data, n.cores, log_file, ...)
   #results <- mergeListDataframes(results)
   return(results)
diff --git a/dev/flat_full.Rmd b/dev/flat_full.Rmd
index 0cd25064ec59cfea270aea87dc56d6531d124c4f..925148be38e0f42345d0a3c660cf4316a1b0c2c4 100644
--- a/dev/flat_full.Rmd
+++ b/dev/flat_full.Rmd
@@ -2437,6 +2437,27 @@ isValidInput2fit <- function(data2fit, formula){
 }
 
 
+
+#' Check if group by exist in data
+#'
+#' @param data The data framecontaining the variables to be used for model fitting.
+#' @param group_by Column name in data representing the grouping variable 
+#'
+#' @return \code{TRUE} if exist otherwise an error is raised
+#'
+#' @examples
+#' is_validGroupBy(mtcars, 'mpg')
+#' @export
+is_validGroupBy <- function(data, group_by){
+  validGroupBy <- group_by %in% names(data)
+  if (!validGroupBy) 
+    stop("<Group by> doen't exist in data !")
+  return(TRUE)
+}
+
+
+
+
 #' Drop Random Effects from a Formula
 #'
 #' This function allows you to remove random effects from a formula by specifying 
@@ -2517,7 +2538,7 @@ is_fullrank <- function(metadata, formula) {
 
 
 #' Fit a model using the fitModel function.
-#' @param group group id to save in glmmTMB obj (usefull for update !)
+#' @param group ID to fit
 #' @param formula Formula specifying the model formula
 #' @param data Data frame containing the data
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
@@ -2525,7 +2546,7 @@ is_fullrank <- function(metadata, formula) {
 #' @export
 #' @examples
 #' fitModel("mtcars" , formula = mpg ~ cyl + disp, data = mtcars)
-fitModel <- function(group, formula, data, ...) {
+fitModel <- function(group , formula, data, ...) {
   # Fit the model using glm.nb from the GLmmTMB package
   model <- glmmTMB::glmmTMB(formula, ..., data = data ) 
   model$frame <- data
@@ -2544,24 +2565,22 @@ fitModel <- function(group, formula, data, ...) {
 
 #' Fit the model based using fitModel functions.
 #'
-#' @param group The specific group to fit the model for
+#' @param groups list of group ID
 #' @param group_by Column name in data representing the grouping variable
-#' @param formula Formula specifying the model formula
 #' @param data Data frame containing the data
-#' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
-#' @return Fitted model object or NULL if there was an error
+#' @return list of dataframe 
 #' @export
+#' @importFrom stats setNames
 #' @examples
-#' subsetData_andfit(group = "setosa", group_by = "Species", 
-#'                  formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
+#' prepare_dataParallel(groups = iris$Species, group_by = "Species", 
 #'                  data = iris )
-subsetData_andfit <- function(group, group_by, formula, data, ...) {
-  subset_data <- data[data[[group_by]] == group, ]
-  fit_res <- fitModel(group, formula, subset_data, ...)
-  #glance_df <- glance.negbin(group_by ,group , fit_res)
-  #tidy_df <- tidy.negbin(group_by ,group,fit_res )
-  #list(glance = glance_df, summary = tidy_df)
-  fit_res
+prepare_dataParallel <- function(groups, group_by, data) {
+  
+  l_data2parallel <- lapply( stats::setNames(groups, groups) , function( group_id ){
+                      subset_data <- data[ data[[ group_by ]] == group_id, ]
+                      return(subset_data)
+                  })
+  return(l_data2parallel)
 }
 
 
@@ -2571,22 +2590,22 @@ subsetData_andfit <- function(group, group_by, formula, data, ...) {
 #' This function fits the model using the specified group, group_by, formula, and data.
 #' It handles warnings and errors during the fitting process and returns the fitted model or NULL if there was an error.
 #'
-#' @param group The specific group to fit the model for
+#' @param data Data frame containing the data
 #' @param group_by Column name in data representing the grouping variable
 #' @param formula Formula specifying the model formula
-#' @param data Data frame containing the data
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
 #' @return List with 'glance' and 'summary' attributes representing the fitted model or NULL if there was an error
 #' @export
 #' @examples
-#' launchFit(group = "setosa", group_by = "Species", 
+#' launchFit(group_by = "Species", 
 #'            formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
-#'            data = iris )
-launchFit <- function(group, group_by, formula, data, ...) {
+#'            data = iris[ iris[["Species"]] == "setosa" , ] )
+launchFit <- function(data, group_by, formula, ...) {
+  group <- unique(data[[ group_by ]]) 
   tryCatch(
     expr = {
       withCallingHandlers(
-          subsetData_andfit(group, group_by, formula, data, ...),
+          fitModel(group , formula, data, ...),
           warning = function(w) {
             message(paste(Sys.time(), "warning for group", group, ":", conditionMessage(w)))
             invokeRestart("muffleWarning")
@@ -2595,7 +2614,6 @@ launchFit <- function(group, group_by, formula, data, ...) {
     error = function(e) {
       message(paste(Sys.time(), "error for group", group, ":", conditionMessage(e)))
       NULL
-      #return(list(glance = empty.glance.negbin(group_by, group), summary = empty.tidy.negbin(group_by, group)))
     }
   )
 }
@@ -2613,10 +2631,9 @@ launchFit <- function(group, group_by, formula, data, ...) {
 #' @param log_file File to write log (default : Rtmpdir/htrfit.log)
 #' @param ... Additional arguments to be passed to the glmmTMB::glmmTMB function
 #' @return List of fitted model objects or NULL for any errors
-#' @importFrom stats setNames
 #' @export
 #' @examples
-#' parallel_fit(group_by = "Species", "setosa", 
+#' parallel_fit(group_by = "Species", groups =  iris$Species, 
 #'                formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
 #'                data = iris, n.cores = 1 )
 parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL, 
@@ -2624,10 +2641,14 @@ parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL,
   
   if (is.null(n.cores)) n.cores <- parallel::detectCores()
   
+  ## get data for parallelization
+  l_data2parallel <- prepare_dataParallel(groups, group_by, data)
+
   clust <- parallel::makeCluster(n.cores, outfile = log_file)
-  parallel::clusterExport(clust, c("subsetData_andfit", "fitModel"),  envir=environment())
-  results_fit <- parallel::parLapply(clust, X = stats::setNames(groups, groups), fun = launchFit, 
-                      group_by = group_by, formula = formula, data = data, ...)
+  parallel::clusterExport(clust, c("fitModel"),  envir=environment())
+  results_fit <- parallel::parLapply(clust, X = l_data2parallel, 
+                                     fun = launchFit, 
+                                     group_by = group_by, formula = formula, ...)
                                      
   parallel::stopCluster(clust)
   #closeAllConnections()
@@ -2651,15 +2672,16 @@ parallel_fit <- function(groups, group_by, formula, data, n.cores = NULL,
 #'                  data = iris, group_by = "Species", n.cores = 1) 
 fitModelParallel <- function(formula, data, group_by, n.cores = NULL, log_file = paste(tempdir(check = FALSE), "htrfit.log", sep = "/"), ...) {
   
-  ## SOme verification
+  ## Some verification
   isValidInput2fit(data, formula)
   is_fullrank(data, formula)
-  
+  is_validGroupBy(data, group_by)
+
   ## -- print log location
   message( paste("Log file location", log_file, sep =': ') ) 
   
-  groups <- unique(data[[group_by]])
   # Fit models in parallel and capture the results
+  groups <- unique(data[[ group_by ]])
   results <- parallel_fit(groups, group_by, formula, data, n.cores, log_file, ...)
   #results <- mergeListDataframes(results)
   return(results)
@@ -2785,83 +2807,34 @@ test_that("Identify full-rank model matrix (with random eff)", {
   expect_true(is_fullrank(metadata, formula))
 })
 
-#test_that(".fitMixteModel returns a fitted mixed-effects model object or NULL if there was an error", {
-#  data(mtcars)
-#  formula <- mpg ~ cyl + disp + (1|gear)
-#  fitted_model <- .fitMixteModel(formula, mtcars)
-  # Add appropriate expectations for the fitted mixed-effects model object
-  
-  # Test with invalid formula
-#  invalid_formula <- formula + "invalid"
-#  fitted_model_error <- .fitMixteModel(invalid_formula, mtcars)
-#  expect_null(fitted_model_error)
-#})
 
-test_that("subsetData_andfit returns a glmTMB obj", {
+test_that("prepare_dataParallel returns a list of dataframe", {
+  ## -- valid input
   data(iris)
-  group <- "setosa"
   group_by <- "Species"
-  formula <- Sepal.Length ~ Sepal.Width + Petal.Length
-  fitted_model <- subsetData_andfit(group, group_by, formula, iris)
-  expect_s3_class(fitted_model, "glmmTMB")
-
-  # Test with invalid formula
-  invalid_formula <- Sepal.Length ~ Sepal.Width + Petal.Length +  invalid_var
-  expect_error(subsetData_andfit(group, group_by, invalid_formula, mtcars))
-  
-  
-    # Additional parameters: 
-   #change family + formula
-  formula <- Sepal.Length ~ Sepal.Width + Petal.Length + (1 | Species)
-  fitted_models <- suppressWarnings(subsetData_andfit(group,
-                                                       group_by,
-                                                       formula = formula, 
-                                                        data = iris, 
-                                                        family = glmmTMB::nbinom1(link = "log") ))
-  expect_s3_class(fitted_models$call$family, "family")
-  expect_equal(fitted_models$call$formula, formula)
-  #change control settings
-  fitted_models <- suppressWarnings(subsetData_andfit(group,
-                                                       group_by,
-                                                       formula = formula, 
-                                                        data = iris, 
-                                                    family = glmmTMB::nbinom1(link = "log"), 
-                                                control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
-                                                                                               eval.max=1e3))))
-  expect_equal(fitted_models$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
+  groups <- unique(iris$Species)
+  l_data <- prepare_dataParallel(groups , group_by, iris)
+  expect_type(l_data, "list")
+  expect_equal(names(l_data), c("setosa", "versicolor", "virginica"))
   
 })
 
 test_that("launchFit handles warnings and errors during the fitting process", {
-  data(mtcars)
-  group <- "Group1"
-  group_by <- "Group"
-  formula <- mpg ~ cyl + disp
-  fitted_model <- suppressWarnings(launchFit(group, group_by, formula, mtcars))
-  expect_s3_class(fitted_model, "glmmTMB")
 
-  # Test with invalid formula
-  invalid_formula <- Sepal.Length ~ Sepal.Width + Petal.Length 
-  output_msg <- capture_message( launchFit(group, group_by, invalid_formula, mtcars))
-  expect_match(output_msg$message, ".* error for group Group1 : object 'Sepal.Length' not found")
-  
-  
   # Additional parameters: 
    #change family + formula
   formula <- Sepal.Length ~ Sepal.Width + Petal.Length
   fitted_models <- suppressWarnings(launchFit(formula = formula, 
-                                                    data = iris, 
-                                                    group_by = group_by, 
-                                                    group = "setosa",
+                                                    data = iris[ iris$Species == "setosa" , ], 
+                                                    group_by = "Species", 
                                                     family = glmmTMB::nbinom1(link = "log") ))
   expect_s3_class(fitted_models$call$family, "family")
   expect_equal(fitted_models$call$formula, formula)
   #change control settings
   fitted_models <- suppressWarnings(launchFit(formula = formula, 
-                                                    data = iris, 
-                                                    group_by = group_by, 
-                                                    group = "setosa",
-                                                     family = glmmTMB::nbinom1(link = "log"), 
+                                                    data = iris[ iris$Species == "setosa" , ], 
+                                                    group_by = "Species", 
+                                                    family = glmmTMB::nbinom1(link = "log") , 
                                                 control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
                                                                                                eval.max=1e3))))
   expect_equal(fitted_models$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
@@ -2911,8 +2884,6 @@ test_that("fitModelParallel fits models in parallel for each group and returns a
   groups <- unique(iris$Species)
   group_by <- "Species"
   formula <- Sepal.Length ~ Sepal.Width + Petal.Length
-  #is.numeric(iris)
-  #iris <- data.frame(lapply(iris, function(y) if(is.numeric(y)) round(y, 0) else y)) 
   fitted_models <- fitModelParallel(formula, iris, group_by, n.cores = 1)
   expect_s3_class(fitted_models$setosa, "glmmTMB")
   expect_length(fitted_models, length(groups))
@@ -2939,6 +2910,12 @@ test_that("fitModelParallel fits models in parallel for each group and returns a
                                                 control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
                                                                                                eval.max=1e3))))
   expect_equal(fitted_models$setosa$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
+  
+  ## -- invalid group by 
+  data(iris)
+  group_by <- "invalid_groupBy"
+  formula <- Sepal.Length ~ Sepal.Width + Petal.Length
+  expect_error(fitModelParallel(formula, iris, group_by, n.cores = 1))
 
 })
 
diff --git a/man/fitModel.Rd b/man/fitModel.Rd
index 56f945e526c9efd1253547e3f3f06047d30acd41..67a588c7d93d429140a79757ecb2d3c1ef50cc59 100644
--- a/man/fitModel.Rd
+++ b/man/fitModel.Rd
@@ -7,7 +7,7 @@
 fitModel(group, formula, data, ...)
 }
 \arguments{
-\item{group}{group id to save in glmmTMB obj (usefull for update !)}
+\item{group}{ID to fit}
 
 \item{formula}{Formula specifying the model formula}
 
diff --git a/man/is_validGroupBy.Rd b/man/is_validGroupBy.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..7af8d6373674d83e209692eadf03ff8f401e97b6
--- /dev/null
+++ b/man/is_validGroupBy.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fitmodel.R
+\name{is_validGroupBy}
+\alias{is_validGroupBy}
+\title{Check if group by exist in data}
+\usage{
+is_validGroupBy(data, group_by)
+}
+\arguments{
+\item{data}{The data framecontaining the variables to be used for model fitting.}
+
+\item{group_by}{Column name in data representing the grouping variable}
+}
+\value{
+\code{TRUE} if exist otherwise an error is raised
+}
+\description{
+Check if group by exist in data
+}
+\examples{
+is_validGroupBy(mtcars, 'mpg')
+}
diff --git a/man/launchFit.Rd b/man/launchFit.Rd
index 4bb44a5cefcbb828a0c49735c7d2b3dbd2128336..92b9ed65f746e91d66f7d94422776738ac50fa58 100644
--- a/man/launchFit.Rd
+++ b/man/launchFit.Rd
@@ -4,17 +4,15 @@
 \alias{launchFit}
 \title{Launch the model fitting process for a specific group.}
 \usage{
-launchFit(group, group_by, formula, data, ...)
+launchFit(data, group_by, formula, ...)
 }
 \arguments{
-\item{group}{The specific group to fit the model for}
+\item{data}{Data frame containing the data}
 
 \item{group_by}{Column name in data representing the grouping variable}
 
 \item{formula}{Formula specifying the model formula}
 
-\item{data}{Data frame containing the data}
-
 \item{...}{Additional arguments to be passed to the glmmTMB::glmmTMB function}
 }
 \value{
@@ -25,7 +23,7 @@ This function fits the model using the specified group, group_by, formula, and d
 It handles warnings and errors during the fitting process and returns the fitted model or NULL if there was an error.
 }
 \examples{
-launchFit(group = "setosa", group_by = "Species", 
+launchFit(group_by = "Species", 
            formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
-           data = iris )
+           data = iris[ iris[["Species"]] == "setosa" , ] )
 }
diff --git a/man/parallel_fit.Rd b/man/parallel_fit.Rd
index 77ebb71c1ed0cb50b1c270953b624558f66ffb43..60fb10ef1dc021361d5b02f6c512a8e226d36a69 100644
--- a/man/parallel_fit.Rd
+++ b/man/parallel_fit.Rd
@@ -39,7 +39,7 @@ Fit models in parallel for each group using mclapply and handle logging.
 Uses parallel_fit to fit the models.
 }
 \examples{
-parallel_fit(group_by = "Species", "setosa", 
+parallel_fit(group_by = "Species", groups =  iris$Species, 
                formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
                data = iris, n.cores = 1 )
 }
diff --git a/man/prepare_dataParallel.Rd b/man/prepare_dataParallel.Rd
new file mode 100644
index 0000000000000000000000000000000000000000..e21e637a0c0c6b4084e087b97d89051608794b92
--- /dev/null
+++ b/man/prepare_dataParallel.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/fitmodel.R
+\name{prepare_dataParallel}
+\alias{prepare_dataParallel}
+\title{Fit the model based using fitModel functions.}
+\usage{
+prepare_dataParallel(groups, group_by, data)
+}
+\arguments{
+\item{groups}{list of group ID}
+
+\item{group_by}{Column name in data representing the grouping variable}
+
+\item{data}{Data frame containing the data}
+}
+\value{
+list of dataframe
+}
+\description{
+Fit the model based using fitModel functions.
+}
+\examples{
+prepare_dataParallel(groups = iris$Species, group_by = "Species", 
+                 data = iris )
+}
diff --git a/man/subsetData_andfit.Rd b/man/subsetData_andfit.Rd
deleted file mode 100644
index e84b1baf65417fc4723a2ebb6f294ebcb6dc1529..0000000000000000000000000000000000000000
--- a/man/subsetData_andfit.Rd
+++ /dev/null
@@ -1,30 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/fitmodel.R
-\name{subsetData_andfit}
-\alias{subsetData_andfit}
-\title{Fit the model based using fitModel functions.}
-\usage{
-subsetData_andfit(group, group_by, formula, data, ...)
-}
-\arguments{
-\item{group}{The specific group to fit the model for}
-
-\item{group_by}{Column name in data representing the grouping variable}
-
-\item{formula}{Formula specifying the model formula}
-
-\item{data}{Data frame containing the data}
-
-\item{...}{Additional arguments to be passed to the glmmTMB::glmmTMB function}
-}
-\value{
-Fitted model object or NULL if there was an error
-}
-\description{
-Fit the model based using fitModel functions.
-}
-\examples{
-subsetData_andfit(group = "setosa", group_by = "Species", 
-                 formula = Sepal.Length ~ Sepal.Width + Petal.Length, 
-                 data = iris )
-}
diff --git a/tests/testthat/test-fitmodel.R b/tests/testthat/test-fitmodel.R
index edd210ae2288e3e291110e3447afe8dc95d1eb45..e6a11fc56206bbd54735b905ccc9369a0b2d3cf4 100644
--- a/tests/testthat/test-fitmodel.R
+++ b/tests/testthat/test-fitmodel.R
@@ -115,83 +115,34 @@ test_that("Identify full-rank model matrix (with random eff)", {
   expect_true(is_fullrank(metadata, formula))
 })
 
-#test_that(".fitMixteModel returns a fitted mixed-effects model object or NULL if there was an error", {
-#  data(mtcars)
-#  formula <- mpg ~ cyl + disp + (1|gear)
-#  fitted_model <- .fitMixteModel(formula, mtcars)
-  # Add appropriate expectations for the fitted mixed-effects model object
-  
-  # Test with invalid formula
-#  invalid_formula <- formula + "invalid"
-#  fitted_model_error <- .fitMixteModel(invalid_formula, mtcars)
-#  expect_null(fitted_model_error)
-#})
 
-test_that("subsetData_andfit returns a glmTMB obj", {
+test_that("prepare_dataParallel returns a list of dataframe", {
+  ## -- valid input
   data(iris)
-  group <- "setosa"
   group_by <- "Species"
-  formula <- Sepal.Length ~ Sepal.Width + Petal.Length
-  fitted_model <- subsetData_andfit(group, group_by, formula, iris)
-  expect_s3_class(fitted_model, "glmmTMB")
-
-  # Test with invalid formula
-  invalid_formula <- Sepal.Length ~ Sepal.Width + Petal.Length +  invalid_var
-  expect_error(subsetData_andfit(group, group_by, invalid_formula, mtcars))
-  
-  
-    # Additional parameters: 
-   #change family + formula
-  formula <- Sepal.Length ~ Sepal.Width + Petal.Length + (1 | Species)
-  fitted_models <- suppressWarnings(subsetData_andfit(group,
-                                                       group_by,
-                                                       formula = formula, 
-                                                        data = iris, 
-                                                        family = glmmTMB::nbinom1(link = "log") ))
-  expect_s3_class(fitted_models$call$family, "family")
-  expect_equal(fitted_models$call$formula, formula)
-  #change control settings
-  fitted_models <- suppressWarnings(subsetData_andfit(group,
-                                                       group_by,
-                                                       formula = formula, 
-                                                        data = iris, 
-                                                    family = glmmTMB::nbinom1(link = "log"), 
-                                                control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
-                                                                                               eval.max=1e3))))
-  expect_equal(fitted_models$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
+  groups <- unique(iris$Species)
+  l_data <- prepare_dataParallel(groups , group_by, iris)
+  expect_type(l_data, "list")
+  expect_equal(names(l_data), c("setosa", "versicolor", "virginica"))
   
 })
 
 test_that("launchFit handles warnings and errors during the fitting process", {
-  data(mtcars)
-  group <- "Group1"
-  group_by <- "Group"
-  formula <- mpg ~ cyl + disp
-  fitted_model <- suppressWarnings(launchFit(group, group_by, formula, mtcars))
-  expect_s3_class(fitted_model, "glmmTMB")
 
-  # Test with invalid formula
-  invalid_formula <- Sepal.Length ~ Sepal.Width + Petal.Length 
-  output_msg <- capture_message( launchFit(group, group_by, invalid_formula, mtcars))
-  expect_match(output_msg$message, ".* error for group Group1 : object 'Sepal.Length' not found")
-  
-  
   # Additional parameters: 
    #change family + formula
   formula <- Sepal.Length ~ Sepal.Width + Petal.Length
   fitted_models <- suppressWarnings(launchFit(formula = formula, 
-                                                    data = iris, 
-                                                    group_by = group_by, 
-                                                    group = "setosa",
+                                                    data = iris[ iris$Species == "setosa" , ], 
+                                                    group_by = "Species", 
                                                     family = glmmTMB::nbinom1(link = "log") ))
   expect_s3_class(fitted_models$call$family, "family")
   expect_equal(fitted_models$call$formula, formula)
   #change control settings
   fitted_models <- suppressWarnings(launchFit(formula = formula, 
-                                                    data = iris, 
-                                                    group_by = group_by, 
-                                                    group = "setosa",
-                                                     family = glmmTMB::nbinom1(link = "log"), 
+                                                    data = iris[ iris$Species == "setosa" , ], 
+                                                    group_by = "Species", 
+                                                    family = glmmTMB::nbinom1(link = "log") , 
                                                 control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
                                                                                                eval.max=1e3))))
   expect_equal(fitted_models$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
@@ -241,8 +192,6 @@ test_that("fitModelParallel fits models in parallel for each group and returns a
   groups <- unique(iris$Species)
   group_by <- "Species"
   formula <- Sepal.Length ~ Sepal.Width + Petal.Length
-  #is.numeric(iris)
-  #iris <- data.frame(lapply(iris, function(y) if(is.numeric(y)) round(y, 0) else y)) 
   fitted_models <- fitModelParallel(formula, iris, group_by, n.cores = 1)
   expect_s3_class(fitted_models$setosa, "glmmTMB")
   expect_length(fitted_models, length(groups))
@@ -269,6 +218,12 @@ test_that("fitModelParallel fits models in parallel for each group and returns a
                                                 control = glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,
                                                                                                eval.max=1e3))))
   expect_equal(fitted_models$setosa$call$control,  glmmTMB::glmmTMBControl(optCtrl=list(iter.max=1e3,eval.max=1e3)))
+  
+  ## -- invalid group by 
+  data(iris)
+  group_by <- "invalid_groupBy"
+  formula <- Sepal.Length ~ Sepal.Width + Petal.Length
+  expect_error(fitModelParallel(formula, iris, group_by, n.cores = 1))
 
 })
 
diff --git a/vignettes/htrfit.Rmd b/vignettes/htrfit.Rmd
index cbd246aac8a390f7da142a7cd89d2d7ff57473b2..f299228a706076ce8b1ef95d4340e96c45e5ae14 100644
--- a/vignettes/htrfit.Rmd
+++ b/vignettes/htrfit.Rmd
@@ -15,7 +15,6 @@ knitr::opts_chunk$set(
 ```
 
 ```{r setup}
-devtools::load_all()
 library(HTRfit)
 ```