R/simulxR.R

In RsSimulx: Extension of 'lixoftConnectors' for 'Simulx'

#' Simulation of mixed effects models and longitudinal data
#'
#' Compute predictions and sample data from \code{Mlxtran} and \code{R} models
#'
#' simulx takes advantage of the modularity of hierarchical models for simulating
#' different components of a model: models for population parameters, individual
#' covariates, individual parameters and longitudinal data.
#'
#' Furthermore, \code{simulx} allows to draw different types of longitudinal data,
#' including continuous, count, categorical, and time-to-event data.
#'
#' The models are encoded using either the model coding language \samp{Mlxtran}.
#' \samp{Mlxtran} models are automatically converted into C++ codes,
#' compiled on the fly and linked to R using the \samp{RJSONIO} package.
#' That allows one to implement very easily complex models and to take advantage
#' of the numerical sovers used by the C++ \samp{mlxLibrary}.
#'
#' See http://simulx.lixoft.com for more details.
#' @param model a \code{Mlxtran} model used for the simulation. It can be a text file or an outut of the inLine function.
#' @param parameter One of
#' \itemize{
#'   \item a vector of parameters with their names and values,
#'   \item a dataframe with parameters defined for each id or each pop
#'   \item a string, path to a data frame (csv or txt file)
#'   \item a string corresponding to the parameter elements automatically generated by monolix
#'   (only when simulation is based on a Monolix project).
#'   One of the following mlx parameter elements: "mlx_Pop", "mlx_PopUncertainSA",
#'   "mlx_PopUncertainLin", "mlx_PopIndiv", "mlx_PopIndivCov","mlx_CondMean",
#'   "mlx_EBEs","mlx_CondDistSample"
#'  }
#' @param covariate One of
#' \itemize{
#'   \item a vector of covariates with their names and values.
#'   \item a dataframe with covariates defined for each id
#'   \item a string, path to a data frame (csv or txt file)
#'   \item a string corresponding to the covariate elements automatically generated by monolix
#'   (only when simulation is based on a Monolix project).
#'   One of the following mlx covariate elements: "mlx_Cov" and "mlx_CovDist"
#' }
#' @param output output or list of outputs. An output can be defined by
#' \itemize{
#'   \item a string corresponding to the output elements automatically generated by monolix
#'   (only when simulation is based on a Monolix project) - the format is mlx_nameofoutput.
#'   \item a list with fields
#'     \itemize{
#'       \item \code{name}: a vector of output names
#'       \item \code{time}:
#'         \itemize{
#'         \item a vector of times
#'         \item a dataframe with columns id, time (columns lloq, uloq, limit are optional)
#'         \item a string, path to a data frame (csv or txt file)
#'       }
#'     \item \code{lloq}: lower limit of quantification (when time is a vector of times)
#'     \item \code{uloq}: upper limit of quantification (when time is a vector of times)
#'     \item \code{limit}: lower bound of the censoring interval (when time is a vector of times)
#'    }
#' }
#' @param treatment treatment or list of treatments. A treatment can be defined by
#' \itemize{
#'   \item A list with fields
#'     \itemize{
#'     \item \code{time} : a vector of input times,
#'     \item \code{amount} : a scalar or a vector of amounts,
#'     \item \code{rate} : a scalar or a vector of infusion rates (default=\code{Inf}),
#'     \item \code{tinf} : a scalar or a vector of infusion times (default=0),
#'     \item \code{washout} : a scalar or a vector of boolean (default=F),
#'     \item \code{adm} : (or type) the administration type (default=1),
#'     \item \code{repeats} : the treatment cycle (optional), a vector with fields
#'       \itemize{
#'       \item \code{cycleDuration} : the duration of a cycle,
#'       \item \code{NumberOfRepetitions} : the number of cycle repetition
#'       }
#'     \item \code{probaMissDose}: the probability to miss each dose (optional).
#'     }
#'  \item a dataframe with treatments defined for each id
#'  \item a string, path to a data frame (csv or txt file)
#'  \item a string corresponding to the treatment elements automatically generated by monolix
#'   (only when simulation is based on a Monolix project) - for example mlx_Adm1.
#' }
#' @param regressor treatment or list of treatments. A treatment can be defined by
#' \itemize{
#'   \item A list with fields
#'     \itemize{
#'     \item \code{name} : a vector of regressor names,
#'     \item \code{time} : a vector of times,
#'     \item \code{value} : a vector of values.
#'     }
#'  \item A dataframe with columns id, time, and regressor values
#'  \item a string, path to a data frame (csv or txt file)
#' }
#' @param occasion An occasion can be defined by
#' \itemize{
#'   \item A list with fields
#'     \itemize{
#'     \item \code{name} : name of the variable which defines the occasions,
#'     \item \code{time} : a vector of times (beginnings of occasions)
#'     }
#'   \item A dataframe with columns id, time, occasion
#'   \item a string, path to a data frame (csv or txt file)
#'   \item "none", to delete an occasion structure
#' }
#' @param varlevel deprecated, use occasion instead.
#' @param group a list, or a list of lists, with fields:
#' \itemize{
#'   \item \code{size} : size of the group (default=1),
#'   \item \code{parameter} : if different parameters per group are defined,
#'   \item \code{covariate} : if different covariates per group are defined,
#'   \item \code{output} : if different outputs per group are defined,
#'   \item \code{treatment} : if different treatments per group are defined,
#'   \item \code{regressor} : if different regression variables per group are defined.
#' }
#' "level" field is not supported anymore in RsSimulx.
#' @param addlines a list with fields:
#' \itemize{
#'   \item \code{formula}: string, or vector of strings, to be inserted .
#' }
#' "section", "block" field are not supported anymore in RsSimulx.
#' You only need to specify a formula. The additional lines will be added in a new section EQUATION.
#' @param project the name of a Monolix project
#' @param nrep Samples with or without uncertainty depending which element is given as "parameter".
#' @param npop deprecated, Set parameter = "mlx_popUncertainSA" or "mlx_popUncertainLin" instead.
#' @param fim deprecated, Set parameter = "mlx_popUncertainSA" or "mlx_popUncertainLin" instead.
#' @param result.file deprecated
#' @param saveSmlxProject If specified, smlx project will be save in the path location
#' (by default smlx project is not saved)
#' @param settings a list of optional settings
#' \itemize{
#'   \item \code{seed} : initialization of the random number generator (integer) (by default a random seed will be generated)
#'   \item \code{id.out}  : add (TRUE) / remove (FALSE) columns id and group when only one element (N = 1 or group = 1) (default=FALSE)
#'   \item \code{kw.max} : deprecated.
#'   \item \code{replacement} : deprecated, use samplingMethod instead
#'   \item \code{samplingMethod}: str, Sampling method used for the simulation. One of "keepOrder", "withReplacement", "withoutReplacement"
#'   (default "keepOrder")
#'   \item \code{out.trt} : TRUE/FALSE (default = TRUE) output of simulx includes treatment
#'   \item \code{out.reg} : TRUE/FALSE (default = TRUE) output of simulx includes regressors
#'   \item \code{sharedIds}: Vector of Elements that share ids. Available types are "covariate", "output",
#'   "treatment", "regressor", "population", "individual" (default c())
#'   \item \code{sameIndividualsAmongGroups}: boolean, if True same individuals will be simulated among all groups
#'   (default False)
#'   \item \code{exportData}: boolean, if True and if a path to save the smlx project (saveSmlxProject) is specified,
#'   export the simulated dataset (smlx project directory/Simulation/simulatedData.txt)
#'   (default False)
#' }
#'
#' @return A list of data frames. Each data frame is an output of simulx
#'
#' @examples
#' \dontrun{
#' myModel <- inlineModel("
#' [LONGITUDINAL]
#' input = {A, k, c, a}
#' EQUATION:
#' t0    = 0
#' f_0   = A
#' ddt_f = -k*f/(c+f)
#' DEFINITION:
#' y = {distribution=normal, prediction=f, sd=a}
#' [INDIVIDUAL]
#' input = {k_pop, omega}
#' DEFINITION:
#' k = {distribution=lognormal, prediction=k_pop, sd=omega}
#' ")
#' 
#' f <- list(name='f', time=seq(0, 30, by=0.1))
#' y <- list(name='y', time=seq(0, 30, by=2))
#' parameter <- c(A=100, k_pop=6, omega=0.3, c=10, a=2)
#' 
#' res <- simulx(model     = myModel,
#'               parameter = parameter,
#'               occasion  = data.frame(time=c(0, 0), occ=c(1, 2)),
#'               output    = list(f,y),
#'               group     = list(size=4),
#'               saveSmlxProject = "./project.smlx")
#' 
#' res <- simulx(model     = myModel,
#'               parameter = parameter,
#'               occasion  = data.frame(time = c(0, 0, 0, 0),
#'                                      occ1 = c(1, 1, 2, 2),
#'                                      occ2 = c(1, 2, 3, 4)),
#'               output    = list(f,y),
#'               group     = list(size=4))
#'               
#' res <- simulx(model     = myModel,
#'               parameter = parameter,
#'               output    = list(f,y),
#'               group     = list(size=4))
#'
#' plot(ggplotmlx() + geom_line(data=res$f, aes(x=time, y=f, colour=id)) +
#'      geom_point(data=res$y, aes(x=time, y=y, colour=id)))
#' print(res$parameter)
#'
#' }
#'
#' @importFrom stats runif
#' @importFrom utils read.table
#' @importFrom gridExtra grid.arrange
#' @export
simulx <- function(model=NULL, parameter=NULL, covariate=NULL, output=NULL,
                   treatment=NULL, regressor=NULL, occasion=NULL, varlevel=NULL,
                   group=NULL, project=NULL, nrep=1, npop=NULL, fim=NULL,
                   saveSmlxProject=NULL, result.file=NULL, addlines=NULL, settings=NULL) {

  params <- as.list(match.call(expand.dots=T))[-1]

  # connectors
  if (!initRsSimulx()$status)
    return()

  # Deprecated Warnings --------------------------------------------------------

  # Message for deprecated varlevel
  if (is.element("varlevel", names(params))) {
    message("[INFO] 'varlevel' argument is deprecated; use 'occasion' instead.")
    if (is.element("occasion", names(params))) {
      stop("varlevel has been replaced by occasion argument. They cannot be used at the same time.", call.=F)
    }
    occasion <- varlevel
  }

  # Message for deprecated result.file
  if (is.element("result.file", names(params))) {
    message("[INFO] 'result.file' argument is deprecated, it will be ignored. ",
            "To save the results, define a path to save the simulx project ",
            "with saveSmlxProject and set the setting exportData to True.")
  }
  
  # Message for deprecated npop
  if (any(c("npop", "fim") %in% names(params))) {
    message("[INFO] 'npop' and 'fim' arguments are deprecated.",
            " You can now define uncertainty with parameters 'mlx_popUncertainSA' or 'mlx_popUncertainLin'",
            " and the number of samples using nrep.")
    # add warning when both npop and parameter are used that npop will be understood as nrep
    if (all(c("parameter", "npop") %in% names(params))) {
      warning("When both parameter and npop are defined, npop will be understood as nrep.", call.=F)
    }
  }

  #=============================================================================
  # Check the inputs of Simulx
  #=============================================================================

  if (!is.null(nrep)) {
    .check_strict_pos_integer(nrep, "nrep")
  }

  # RETRO 2020
  if (!is.null(npop)) {
    .check_strict_pos_integer(npop, "npop")
    runSimPop <- TRUE
  } else {
    runSimPop <- FALSE
  }

  # Model file -----------------------------------------------------------------
  nbID <- 1
  if(!is.null(model)){
    .checkModel(fileName = model)
  }
  # set_options(warnings = FALSE)

  # Settings -------------------------------------------------------------------
  if (is.element("kw.max", names(settings))) {
    if (is.null(npop) | ! is.null(parameter)) {
      message("[INFO] Setting kw.max is deprecated, it will be ignored (it can only be used with simpopmlx).")
    }
  }
  .check_in_vector(
    names(settings), "settings",
    c("seed", "kw.max", "sep", "digits", "replacement", "samplingMethod", "out.trt",
      "out.reg", "id.out", "sameIndividualsAmongGroups", "sharedIds", "exportData")
  )
  settings <- .initsimulxSettings(settings)

  # Output files ---------------------------------------------------------------
  if (!is.null(saveSmlxProject)) {
    .check_char(saveSmlxProject, "saveSmlxProject")
    if (is.null(.getFileExt(saveSmlxProject))) {
      saveSmlxProject <- paste0(saveSmlxProject, ".smlx")
    }
    ext <- .getFileExt(saveSmlxProject)
    if (ext != "smlx") {
      stop("Invalid saveSmlxProject '", saveSmlxProject, "'. File extension must be smlx.", call.=F)
    }
    if (! dir.exists(dirname(saveSmlxProject))) {
      dir.create(dirname(saveSmlxProject), recursive = TRUE)
    }
  }
  
  if (settings$exportData && is.null(saveSmlxProject)) {
    warning("You first need to specify a path to save smlx project in order to export data.", call.=F)
    settings$exportData <- F
  }
  # Project file ---------------------------------------------------------------
  if(!is.null(project)){
    .check_file(filename = project, fileType = "Project")
    if (!.getFileExt(project) == "mlxtran") {
      stop("Invalid project file. Monolix project must have a .mlxtran extension.", call. = FALSE)
    }
  }

  # Project file and model file ------------------------------------------------
  if(is.null(model) & is.null(project)) {
    stop("You must define either the model or the Monolix project file.", call. = FALSE)
  }
  if (!is.null(model) & !is.null(project)) {
    stop("You must define either the model or the Monolix project file, not both of them.", call. = FALSE)
  }

  mlxProject <- !is.null(project)

  # Parameter ------------------------------------------------------------------
  # merge parameters
  parameter <- .transformParameter(parameter)
  # check parameter format
  parameter <- .checkParameter(parameter, mlxProject=mlxProject)

  # Write the data set as a .txt file
  if (is.data.frame(parameter)) {
    parameter <- .addDataFrameTemp(df = parameter)
  }
  # update the size of the default group
  if (is.string(parameter) && file.exists(parameter)) {
    df <- utils::read.table(file=parameter, header=T, sep=.getDelimiter(parameter))
    nbID <- max(nbID, .getNbIds(df))

    # RETRO 2020
    indexPop <- which(names(df) == 'pop')
    if (length(indexPop) > 0) {
      npop <- length(unique(df[,indexPop[1]]))
      runSimPop <- FALSE
    }
  }
  
  # covariates -----------------------------------------------------------------
  # merge covariates
  covariate <- .transformParameter(covariate)
  # check covariate format
  covariate <- .checkCovariate(covariate, mlxProject=mlxProject)
  # Write the data set as a .txt file
  if (is.data.frame(covariate)) {
    covariate <- .addDataFrameTemp(df = covariate)
  }
  # update the size of the default group
  if (is.string(covariate) && file.exists(covariate)) {
    df <- utils::read.table(file=covariate, header=T, sep=.getDelimiter(covariate))
    nbID <- max(nbID, .getNbIds(df))
  }
    
  # Output ---------------------------------------------------------------------
  # check output format
  output <- .checkOutput(output, mlxProject=mlxProject)
  # if dataframe, save it in a txt file
  for (iout in seq_along(output)) {
    if ("time" %in% names(output[[iout]]) && is.data.frame(output[[iout]]$time)) {
      outputValue <- output[[iout]]$time
      # Write the data set as a .txt file
      output[[iout]]$time <- .addDataFrameTemp(df=outputValue)
      # update the size of the default group
      nbID <- max(nbID, .getNbIds(outputValue))
    }
  }

  # Treatment ------------------------------------------------------------------
  # check treatment format
  treatment <- .checkTreatment(treatment, mlxProject=mlxProject)
  treatment <- .splitTreatment(treatment)

  for (itreat in seq_along(treatment)) {
    treatementValue <- treatment[[itreat]]
    if (is.data.frame(treatementValue)) {
      # Write the data set as a .txt file
      treatment[[itreat]] <-.addDataFrameTemp(df = treatementValue)
      # update the size of the default group
      nbID <- max(nbID, .getNbIds(treatementValue))
    }
  }

  # Regressor ------------------------------------------------------------------
  # check regressor format
  regressor <- .checkRegressor(regressor)
  # merge regressors
  regressor <- .transformRegressor(regressor)
  # When regressor is defined as dataframe, write the data set as a .txt file
  if (is.data.frame(regressor)) {
    nbID <- max(nbID, .getNbIds(regressor))
    regressor <- .addDataFrameTemp(df = regressor)
  }
  
  # Occasions ------------------------------------------------------------------
  # check occasions format
  occasion <- .checkOccasion(occasion)

  if (! is.null(occasion) && (any(c("id", "ID") %in% names(occasion)))) {
    # Write the data set as a .txt file
    # occasion <- .renameColumns(occasion, "id", "ID")
    occasion <-.addDataFrameTemp(df = occasion)
    # update the size of the default group
    nbID <- max(nbID, .getNbIds(occasion))
  }
  
  # Group ----------------------------------------------------------------------
  allowedNames <- c("size", "parameter", "covariate", "output", "treatment", "regressor", "level")
  if (any(is.element(allowedNames, names(group))) | !all(sapply(group, .is_list_or_named_vector))) {
    group <- list(group)
  }
  groupsWithLevel <- group[sapply(group, function(g) is.element("level", names(g)))]
  for (g in groupsWithLevel) {
    if (any(g[["level"]] != "individual")) {
      stop("'level' information in 'group' argument is not supported by RsSimulx. \n",
           "Simulations are only run at the individual level. If you have variability at other levels, check the online documentation.", call. = FALSE)
    }
  }

  defaultGroup <- list()
  for (element in c("parameter", "covariate", "output", "treatment", "regressor")) {
    defaultGroup[[element]] <- element %in% names(params)
  }
  group <- .checkGroup(group, defaultGroup, mlxProject=mlxProject)

  #=============================================================================
  # Initialize the simulation
  #=============================================================================
  if (!is.null(model)) {
    # The simulation is initialized by a model

    # Check that the output is not null
    isOutputArg <- ! is.null(output) || any(sapply(group, function(g) "output" %in% names(g)))
    isOutputModel <- .checkModelOutputSection(model)

    if (!isOutputArg && !isOutputModel) {
      stop("When a model is defined without an OUTPUT section, you must define an output element ",
           "(either in output or in group definition).", call. = F)
    }
    if (!isOutputModel) {
      outputName <- .getOutputNames(output, group)
      model <- .addOutputToFile(model, outputName[1])
    }

    # create project with the model file
    .lixoftCall("newProject", list(modelFile = model))

  } else {
    # The simulation is initialized by a Monolix project

    # RETRO 2020 - population parameters with simpopmlx
    if (runSimPop) {
      popParam <- simpopmlx(n = npop, project = project, fim = fim,
                            kw.max = settings[["kw.max"]], seed = settings[["seed"]])
      if (is.element("pop", names(popParam))) {
        popParam <- popParam[names(popParam) != "pop"]
      }
    }

    # create project with monolix project
    .lixoftCall("importProject", list(projectFile = project))

    sharedGroupName <- .lixoftCall("getGroups")[[1]]$name

    if (is.element("mlx_PopInit", names(.lixoftCall("getPopulationElements")))) {
      stop("Population parameter estimation has not been launched in Monolix project.", call. = FALSE)
    }
    
    # RETRO 2020 - define population parameters with simpopmlx
    if (runSimPop) {
      .lixoftCall("definePopulationElement", list(name = "popSimpopMlx", element = .addDataFrameTemp(df = popParam)))
      .lixoftCall("setGroupElement", list(group=sharedGroupName, elements = "popSimpopMlx"))
    }

  }

  sharedGroupName <- .lixoftCall("getGroups")[[1]]$name

  # RETRO 2020
  if (!is.null(npop)) {
    nrep <- npop * nrep
  }

  # Occasions
  if (!is.null(occasion)) {
    if (is.string(occasion) && occasion == "none") {
      .lixoftCall("deleteOccasionElement")
    } else {
      .lixoftCall("defineOccasionElement", list(element = occasion))
    }
  }

  # set groups settings size ---------------------------------------------------
  # shared group size
  if (nbID > 1) {
    .lixoftCall("setGroupSize", list(group = sharedGroupName, size = nbID))
  }

  # Add additional lines in the model ------------------------------------------
  if (is.element("formula", names(addlines))) {
    addlines <- list(addlines)
  }
  if (!is.null(addlines)) {
    for (newline in addlines) {
      .checkAdditionalLine(newline)
      .lixoftCall("setAddLines", list(lines = newline$formula))
    }
  }

  #=============================================================================
  # Design the Shared group
  #=============================================================================
  inputParameter <- parameter

  # Add output -----------------------------------------------------------------
  .addOutput(output)

  # Add parameter --------------------------------------------------------------
  if (!is.null(parameter)) {
    if (is.string(parameter) && parameter %in% .getAllowedMlxParameterNames()) {
      .addMlxParameter(parameter)
    } else {
      .addParameter(parameter)
    }
  }

  # Add covariate --------------------------------------------------------------
  # RETRO 2020 - Include covariates from parameter argument
  if (! is.null(parameter) && ! parameter[1] %in% .getAllowedMlxParameterNames()) {
    covFromParameter <- .filterParameter(parameter, "cov", store_dataframe=F, unlistOutput=T)
    if (! is.null(covFromParameter)) {
      message("[INFO] Using 'parameter' argument to define covariates is deprecated. Use 'covariate' argument instead.")
      if (is.string(covariate)) {
        if (covariate %in% .getAllowedMlxCovariateNames()) {
          stop("You must either define covariate with a string corresponding to the element generated by Monolix, or with a list.", call.=F)
        }
        covariate <- utils::read.table(file=covariate, header=T, sep=.getDelimiter(covariate))
      }
      covariate <- .transformParameter(list(
        covFromParameter,
        covariate
      ))
      
      if (is.data.frame(covariate)) {
        covariate <- .addDataFrameTemp(df = covariate)
      }
    }
    if (is.null(c(.filterParameter(parameter, "ind"), .filterParameter(parameter, "pop")))) {
      parameter <- NULL
    }
  }

  if (! is.null(covariate)) {
    if (covariate[1] %in% .getAllowedMlxCovariateNames()) {
      .addMlxCovariate(covariate)
    } else {
      .addCovariate(covariate)
    }
  }

  # Add treatment --------------------------------------------------------------
  if (!is.null(treatment)) {
    .addTreatment(treatment)
  }

  # Add regressor --------------------------------------------------------------
  .addRegressor(regressor)

  groupOut <- NULL

  #=============================================================================
  # Design the groups
  #=============================================================================
  if (!is.null(group)) {
    ngroup <- length(group)
    groupNames <- paste0("simulationGroup", seq(ngroup))
    groupNames[1] <- sharedGroupName

    for (igroup in seq_along(group)) {
      gname <- groupNames[igroup]
      g <- group[[igroup]]

      # group not yet initialized
      if (! gname %in% sapply(.lixoftCall("getGroups"), function(g) g$name)) {
        .lixoftCall("addGroup", list(group=gname))
      }
      
      # fill group information
      group[[igroup]] <- .addGroup(group=g, groupName=gname)

      groupOut[[igroup]] <- list(size=.lixoftCall("getGroups")[[igroup]]$size, level='longitudinal')
    }
    
    if (ngroup == 1) groupOut <- NULL
  }

  # Check missing elements in definition ---------------------------------------
  # Check if output have been defined
  if (! mlxProject && is.null(output) && ! any(sapply(group, function(g) "output" %in% names(g)))) {
    o <- .lixoftCall("getGroups")[[1]]$output
    warning("Output has not been defined. Default output element '", o,
            "' will be used.", call.=F)
  }

  # Check if parameter have been defined
  if (! mlxProject && is.null(parameter) && ! any(sapply(group, function(g) "parameter" %in% names(g)))) {
    p <- .lixoftCall("getGroups")[[1]]$parameter$name
    warning("Parameter has not been defined. Default parameter element '", p,
            "' will be used.", call.=F)
  }
  
  # Check if covariates have been defined
  pType <- .lixoftCall("getGroups")[[1]]$parameter$type
  if (! mlxProject && pType != "individual" && is.null(covariate) && ! any(sapply(group, function(g) "covariate" %in% names(g)))) {
    covElement <- .lixoftCall("getCovariateElements")
    if (length(covElement) >= 1) {
      warning("Covariate has not been defined. Default covariate element '", names(covElement),
              "' will be used.", call.=F)
    }
  }
  
  # Check if regressor have been defined
  if (! mlxProject && is.null(regressor) && ! any(sapply(group, function(g) "regressor" %in% names(g)))) {
    regElement <- .lixoftCall("getRegressorElements")
    if (length(regElement) >= 1) {
      warning("Regressor has not been defined. Default regressor element '", names(regElement),
              "' will be used.", call.=F)
    }
  }
  
  #=============================================================================
  # Manage the replicates
  #=============================================================================
  if (nrep > 1) {
    .lixoftCall("setNbReplicates", list(nb = nrep))
  }

  #=============================================================================
  # Manage the settings
  #=============================================================================
  # shared ids
  if (!is.null(settings$sharedIds)) {
    .lixoftCall("setSharedIds", list(settings$sharedIds))
  }
  
  # sameIndividualsAmongGroups
  if (settings$sameIndividualsAmongGroups == T) {
    if (length(group) <= 1) {
      warning("There is only one shared group, sameIndividualsAmongGroups is set to FALSE", call.=F)
      settings$sameIndividualsAmongGroups <- F
    }
    cov_groups <- sapply(.lixoftCall("getGroups"), function(g) g$covariate)
    if (length(unique(cov_groups)) > 1) {
      warning("Different covariates have been defined between groups, sameIndividualsAmongGroups is set to FALSE", call.=F)
      settings$sameIndividualsAmongGroups <- F
    }
  }
  .lixoftCall("setSameIndividualsAmongGroups", list(settings$sameIndividualsAmongGroups))

  # seed
  .lixoftCall("setProjectSettings", list(seed = settings[['seed']]))
  
  # sampling method
  .lixoftCall("setSamplingMethod", list(method = settings$samplingMethod))

  # rename groups
  for (index in seq_along(.lixoftCall("getGroups"))) {
    .lixoftCall("renameGroup", list(paste0('simulationGroup', index), paste0(index)))
  }

  #=============================================================================
  # Run the simulation
  #=============================================================================
  .lixoftCall("runSimulation")

  simOut <- .lixoftCall("getSimulationResults")$res

  # Manage Simulation Output in results ----------------------------------------
  for (index in seq_along(simOut)) {
    outName <- names(simOut)[index]
    out <- simOut[[index]]

    # censor output
    out <- .censorOutput(out, outName, output, group)

    # factor rep id group cens pop & remove unused columns
    out <- .postProcessDataFrames(out, force=settings[["id.out"]])

    # attributes
    out <- out[c(setdiff(names(out), outName), outName)]
    attributes(out) <- c(attributes(out), name = outName)

    simOut[[index]] <- out
  }
  if(!is.null(groupOut)){
    simOut$group <- groupOut
  }

  # Manage Parameters Output in results ----------------------------------------
  paramOut <- .lixoftCall("getSimulationResults")$IndividualParameters
  if (length(paramOut) == 1) {
    paramOut <- paramOut[[1]]
  } else {
    # groups
    for (g in names(paramOut)) {
      paramOut[[g]] <- cbind(paramOut[[g]], group=g)
    }
    paramOut <- do.call(rbind, paramOut)
    row.names(paramOut) <- seq_len(nrow(paramOut))
  }

  if (ncol(paramOut) > 1) {
    for (indexCol in 2:ncol(paramOut)) {
      notDouble <- suppressWarnings(sum(is.na(as.double(paramOut[,indexCol]) == paramOut[,indexCol])))
      if (notDouble == 0) {
        paramOut[,indexCol] <- as.double(paramOut[,indexCol])
      }
    }
  }

  # factor rep id group cens pop & remove unused columns
  paramOut <- .postProcessDataFrames(paramOut, force = settings[["id.out"]])

  # if only one set of parameter, return parameter as named vector
  if (nrow(unique(paramOut[! names(paramOut) %in% c("rep", "id", "group")])) == 1) {
    paramOut <- utils::head(paramOut[! names(paramOut) %in% c("rep", "id", "group")], n = 1)
  }
  if (nrow(paramOut) == 1) {
    paramOut <- unlist(paramOut)
  }

  simOut$parameter <- paramOut

  # Manage treatment Output in results ----------------------------------------
  if (settings[['out.trt']]) {
    trtFile <- file.path(.lixoftCall("getProjectSettings")$directory, '/Simulation/doses.txt')
    if (file.exists(trtFile)) {
      treat <- read.table(file=trtFile, header=T, sep=.getDelimiter(trtFile))

      # rename column amt to amount
      treat <- .renameColumns(treat, "amt", "amount")

      if (sum(is.na(treat$tinf)) == nrow(treat)) {
        treat <- treat[, names(treat) != "tinf"]
      }
      if (is.element("tinf", names(treat))) {
        treat$tinf[!is.na(treat$tinf)] <- treat$amount[!is.na(treat$tinf)] / treat$tinf[!is.na(treat$tinf)]
        treat <- .renameColumns(treat, "tinf", "rate")
      }

      # factor rep id group cens pop & remove unused columns
      treat <- .postProcessDataFrames(treat, force = settings[["id.out"]])

      simOut$treatment <- treat
    }
  }
  
  # Manage regressor Output in results ----------------------------------------
  if (settings[['out.reg']]) {
    regFile <- file.path(.lixoftCall("getProjectSettings")$directory, '/Simulation/regressors.txt')
    if (file.exists(regFile)) {
      regressors <- read.table(file=regFile, header=T, sep=.getDelimiter(regFile))

      # factor rep id group cens pop & remove unused columns
      regressors <- .postProcessDataFrames(regressors, force = settings[["id.out"]])
      
      simOut$regressors <- regressors
    }
  }

  # Manage population Output in results ----------------------------------------
  if (!is.null(nrep)) {
    popFile <- file.path(.lixoftCall("getProjectSettings")$directory,
                         "Simulation", "populationParameters.txt")

    if (file.exists(popFile)) {
      popData <- read.table(file=popFile, header=T, sep=.getDelimiter(popFile))

      # factor rep id group cens pop & remove unused columns
      popData <- .postProcessDataFrames(popData, force=settings[["id.out"]])

      # if only one set of parameter, return only the first row
      if (nrow(unique(popData[! names(popData) %in% c("rep", "id", "group")])) == 1) {
        popData <- popData[! names(popData) %in% c("rep", "id", "group")][1,]
      }

      # if one row: return a named vector
      if (nrow(popData) == 1) {
        popData <- unlist(popData)
      }

      simOut$population <- popData
    }
  }

  # Write data -----------------------------------------------------------------
  if (!is.null(saveSmlxProject)) {
    .lixoftCall("setProjectSettings", list(list(userfilesnexttoproject=T)))
    .lixoftCall("saveProject", list(saveSmlxProject))
  }

  if (settings$exportData) {
    version <- .lixoftCall("getLixoftConnectorsState")$version
    # RETRO 2020
    if (version == "2020R1") {
      file <- file.path(.lixoftCall("getProjectSettings")$directory, "Simulation", "simulatedData.txt")
      writeData(filename=file)
    } else {
      .lixoftCall("exportSimulatedData")
    }
  }

  return(simOut)
}

.initsimulxSettings <- function(settings) {
  if (is.null(settings)) settings <- list()
  
  # Message for deprecated settings
  if (is.element("replacement", names(settings))) {
    message("[INFO] Setting replacement is deprecated, use samplingMethod setting instead.")
    if (is.element("samplingMethod", names(settings))) {
      stop("replacement has been replaced by samplingMethod setting. They cannot be used at the same time.", call.=F)
    }
    if (settings$replacement == T) {
      settings$samplingMethod <- "withReplacement"
    } else {
      settings$samplingMethod <- "keepOrder"
    }
  }
  if (is.element("digits", names(settings))) {
    message("[INFO] Setting digits is deprecated, it will be ignored.")
  }
  if (is.element("sep", names(settings))) {
    message("[INFO] Setting sep is deprecated, it will be ignored.")
  }
  
  if (!is.element("seed", names(settings))) settings$seed <- sample.int(1e6, 1)
  if (!is.element("kw.max", names(settings))) settings[["kw.max"]] <- 100
  if (!is.element("replacement", names(settings))) settings[['replacement']] <- F
  if (!is.element("samplingMethod", names(settings))) settings[['samplingMethod']] <- "keepOrder"
  if (!is.element("out.trt", names(settings))) settings[['out.trt']] <- T
  if (!is.element("out.reg", names(settings))) settings[['out.reg']] <- T
  if (!is.element("id.out", names(settings))) settings[['id.out']] <- F
  if (!is.element("sharedIds", names(settings))) settings[["sharedIds"]] <- c()
  if (!is.element("sameIndividualsAmongGroups", names(settings))) settings[["sameIndividualsAmongGroups"]] <- F
  if (!is.element("exportData", names(settings))) settings$exportData <- F

  .check_pos_integer(settings$seed, "setting seed")
  .check_pos_integer(settings$kw.max, "setting kw.max")
  .check_bool(settings$replacement, "setting replacement")
  .check_in_vector(settings$samplingMethod, "setting samplingMethod", c("keepOrder", "withReplacement", "withoutReplacement"))
  .check_bool(settings$out.trt, "setting out.trt")
  .check_bool(settings$out.reg, "setting out.reg")
  .check_bool(settings$id.out, "setting id.out")
  .check_in_vector(settings$sharedIds, "setting sharedIds",
                   c("covariate", "output", "treatment", "regressor", "population", "individual"))
  .check_bool(settings$sameIndividualsAmongGroups, "setting sameIndividualsAmongGroups")
  .check_bool(settings$exportData, "setting exportData")

  return(settings)
}

.rearrangeRepPop <- function(df, npop) {
  if (! "rep" %in% names(df)) return(df)
  if (nrow(df) == 0) return(df)
  if (is.null(npop)) npop <- 1
  if (npop == 1) return(df)
  nrep <- length(unique(df$rep)) / npop
  count <- table(df$rep)[[1]]
  pop <- rep(rep(seq_len(npop), each = count), nrep)
  rep <- rep(seq_len(nrep), each = count * npop)
  df$rep <- rep
  df$pop <- pop
  # reorder
  df <- df[, c("rep", "pop", names(df)[! names(df) %in% c("rep", "pop")])]
  return(df)
}

.postProcessDataFrames <- function(df, force = FALSE) {
  cols <- intersect(c("cens", "pop", "rep", "id", "group"), names(df))
  for (col in cols) {
    df[[col]] <- as.factor(df[[col]])
  }

  if (is.element("rep", names(df))) {
    if (length(unique(df$rep)) == 1) df <- df[names(df) != "rep"]
  }

  if (is.element("pop", names(df))) {
    if (length(unique(df$pop)) == 1) df <- df[names(df) != "pop"]
  }

  if (is.element("id", names(df))) {
    if (length(unique(df$id)) == 1 & ! force) df <- df[names(df) != "id"]
  }

  if (is.element("group", names(df))) {
    if (length(unique(df$group)) == 1 & ! force) df <- df[names(df) != "group"]
  }
  return(df)
}