diff --git a/DESCRIPTION b/DESCRIPTION
index d1f1985..7e07728 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: PKPDmap
 Type: Package
 Title: MAP Bayesian estimates
-Version: 1.0.3
-Date: 2025-06-18
+Version: 1.1.4
+Date: 2025-07-12
 Author: Ron Keizer, Jasmine Hughes, Kara Woo
 Maintainer: Ron Keizer <ron@insight-rx.com>
 Description: Obtain MAP Bayesian estimates based on individual data and
diff --git a/R/calc_residuals.R b/R/calc_residuals.R
index 46e3b19..b3c4173 100644
--- a/R/calc_residuals.R
+++ b/R/calc_residuals.R
@@ -18,6 +18,7 @@ calc_residuals <- function(
   parameters_population,
   covariates,
   regimen,
+  lagtime,
   omega_full,
   error,
   weights,
@@ -60,6 +61,7 @@ calc_residuals <- function(
       iov_bins = iov_bins,
       output_include = output_include,
       t_init = t_init,
+      lagtime = lagtime,
       ...
     )
   })
@@ -79,6 +81,7 @@ calc_residuals <- function(
       iov_bins = iov_bins,
       A_init = A_init_population,
       t_init = t_init,
+      lagtime = lagtime,
       ...
     )
   })
diff --git a/R/get_individual_parameters_from_fit.R b/R/get_individual_parameters_from_fit.R
index 25c7a8e..8be8792 100644
--- a/R/get_individual_parameters_from_fit.R
+++ b/R/get_individual_parameters_from_fit.R
@@ -1,7 +1,8 @@
 #' Get the individual parameters from the fitted coefficients (etas)
 #' 
+#' @inheritParams get_map_estimates
 #' @param fit fit object
-#' @param parameters list of population parameters
+#' @param nonfixed vector of parameters that are not fixed (i.e. estimated)
 #'
 get_individual_parameters_from_fit <- function(
   fit,
diff --git a/R/get_map_estimates.R b/R/get_map_estimates.R
index 046e833..696c8ec 100755
--- a/R/get_map_estimates.R
+++ b/R/get_map_estimates.R
@@ -18,6 +18,8 @@
 #' PKPDsim.
 #' @param error residual error, specified as list with arguments `add` and/or 
 #' `prop` specifying the additive and proportional parts
+#' @param lagtime vector of lagtimes for each compartment in the model, either
+#' numeric vector, or vector of parameters.
 #' @param ltbs log-transform both sides? (`NULL` by default, meaning that it 
 #' will be picked up from the PKPDsim model. Can be overridden with `TRUE`). 
 #' Note: `error` should commonly only have additive part.
@@ -39,7 +41,12 @@
 #' @param ll_func likelihood function, default is `ll_func_PKPDsim` as included 
 #' in this package.
 #' @param optimizer optimization library to use, default is `optim`
-#' @param method optimization method, default `BFGS`
+#' @param method optimization method, default `BFGS`. This method is a
+#' gradient-based method, in which the gradients are computed using finite-
+#' difference method. When the model contains a step-function where a
+#' estimated parameter is involved in the step-function (such as in the case of
+#' lagtime), it is advised to use a non-gradient-based estimation method
+#' such as `Nelder-Mead` to avoid estimation failures.
 #' @param control list of options passed to `optim()` function
 #' @param allow_obs_before_dose allow observation before first dose?
 #' @param type estimation type, options are `map`, `ls`, and `np_hybrid`
@@ -112,6 +119,7 @@ get_map_estimates <- function(
                       include_omega = TRUE,
                       include_error = TRUE,
                       regimen = NULL,
+                      lagtime = NULL,
                       t_init = 0,
                       int_step_size = 0.1,
                       ll_func = ll_func_PKPDsim,
@@ -137,6 +145,9 @@ get_map_estimates <- function(
   if(weight_prior_var == 0) {
     calc_ofv <- calc_ofv_ls
   }
+
+  ## Make sure lagtime is properly formatted (potentially pulled from model)
+  lagtime <- PKPDsim:::parse_lagtime(lagtime, model, parameters)
   
   ## Misc checks:
   check_inputs(model, data, parameters, omega, regimen, censoring, type)
@@ -211,6 +222,7 @@ get_map_estimates <- function(
       n_ind = 1,
       int_step_size = int_step_size,
       regimen = regimen,
+      lagtime = lagtime,
       t_obs = t_obs,
       obs_type = data$obs_type,
       checks = FALSE,
@@ -250,6 +262,7 @@ get_map_estimates <- function(
           t_obs = t_obs,
           model = model,
           regimen = regimen,
+          lagtime = lagtime,
           error = error,
           omega_full = omega$est / weight_prior_var,
           omega_inv = solve(omega$est / weight_prior_var),
@@ -305,6 +318,7 @@ get_map_estimates <- function(
           t_obs = t_obs,
           model = model,
           regimen = regimen,
+          lagtime = lagtime,
           error = error,
           nonfixed = omega$nonfixed,
           transf = transf,
@@ -353,6 +367,7 @@ get_map_estimates <- function(
       obj = obj,
       model = model,
       regimen = regimen,
+      lagtime = lagtime,
       data = data,
       covariates = covariates,
       weights = weights,
@@ -399,6 +414,7 @@ get_map_estimates <- function(
       covariates = covariates,
       int_step_size = int_step_size,
       regimen = regimen,
+      lagtime = lagtime,
       omega_full = omega$full,
       error = error,
       weights = weights,
diff --git a/R/ll_func_PKPDsim.R b/R/ll_func_PKPDsim.R
index 8028aa1..cd7f26e 100644
--- a/R/ll_func_PKPDsim.R
+++ b/R/ll_func_PKPDsim.R
@@ -1,29 +1,19 @@
 #' Likelihood function for MAP optimization using PKPDsim
 #' 
-#' @param data vector of
+#' @inheritParams get_map_estimates
 #' @param sim_object design (event-table) obtained from PKPDsim to be used in simulations/optimizations
-#' @param parameters parameter list
-#' @param covariates covariates list
 #' @param nonfixed non-fixed (i.e. estimated) parameters
-#' @param error error model to use, e.g. `list(add = .5, prop = .15)`
-#' @param model PKPDsim model
 #' @param omega_full full omega matrix
 #' @param omega_inv inverse of omega matrix. Passed to this function to avoid doing computations in each iteration of the search.
 #' @param omega_eigen eigenvalue decomposation (logged) of omega matrix. Passed to this function to avoid doing computations in each iteration of the search.
 #' @param sig signficance, used in optimization
-#' @param weights weights for data, generally a vector of weights between 0 and 1.
 #' @param transf transformation function for data, if needed. E.g. `log(x)`. Default is `function(x) = x`.
-#' @param as_eta implement as regular eta instead of exponential eta, can be vector.
 #' @param censoring_idx censoring indices, used for <LOQ data.
 #' @param censoring_label censoring label, used for <LOQ data
-#' @param t_init init time for simulations, default 0.
-#' @param iov_bins IOV bins object
 #' @param calc_ofv function to calculate OFV based on simulated data and set of parameters and omega matrix
-#' @param regimen PKPDsim regimen
 #' @param steady_state_analytic list object with settings for steady state MAP estimation.
 #' @param include_omega include omega in calculation of OFV?
 #' @param include_error include residual error in calculation of OFV?
-#' @param verbose verbose output?
 #' @param eta01 eta1 
 #' @param eta02 eta2 
 #' @param eta03 eta3 
@@ -75,6 +65,7 @@ ll_func_PKPDsim <- function(
   t_init = 0,
   calc_ofv,
   regimen,
+  lagtime = c(0),
   steady_state_analytic,
   include_omega,
   include_error,
@@ -109,7 +100,9 @@ ll_func_PKPDsim <- function(
     sim_object,
     ode = model,
     duplicate_t_obs = TRUE,
-    t_init = t_init)$y)
+    t_init = t_init,
+    lagtime = lagtime
+  )$y)
   dv <- transf(data$y)
   obs_type <- data$obs_type
   ofv_cens <- NULL
diff --git a/R/ll_func_generic.R b/R/ll_func_generic.R
index e739a82..3dce43e 100755
--- a/R/ll_func_generic.R
+++ b/R/ll_func_generic.R
@@ -2,6 +2,34 @@
 #' 
 #' @inheritParams ll_func_PKPDsim
 #'
+#' @param omega_full full omega matrix
+#' @param sig signficance, used in optimization
+#' @param eta01 eta1 
+#' @param eta02 eta2 
+#' @param eta03 eta3 
+#' @param eta04 eta4 
+#' @param eta05 eta5 
+#' @param eta06 eta6 
+#' @param eta07 eta7 
+#' @param eta08 eta8 
+#' @param eta09 eta9 
+#' @param eta10 eta10 
+#' @param eta11 eta11 
+#' @param eta12 eta12 
+#' @param eta13 eta13 
+#' @param eta14 eta14 
+#' @param eta15 eta15 
+#' @param eta16 eta16 
+#' @param eta17 eta17 
+#' @param eta18 eta19
+#' @param eta19 eta19
+#' @param eta20 eta20
+#' @param eta21 eta21
+#' @param eta22 eta22
+#' @param eta23 eta23
+#' @param eta24 eta24 
+#' @param ... passed on to PKPDsim
+#' 
 ll_func_generic <- function(
   data,
   # unfortunately seems no other way to do this...
@@ -11,6 +39,7 @@ ll_func_generic <- function(
   parameters,
   covariates = NULL,
   regimen = regimen,
+  lagtime = NULL,
   omega_full = omega_full,
   error = error,
   model,
diff --git a/R/parse_input_data.R b/R/parse_input_data.R
index 62bd2fd..aefe23b 100644
--- a/R/parse_input_data.R
+++ b/R/parse_input_data.R
@@ -7,6 +7,12 @@ parse_input_data <- function(
   cols = list(x = "t", y = "y"),
   obs_type_label = NULL
 ) {
+  ## Parse data to include label for obs_type
+  if(is.null(obs_type_label)) {
+    data$obs_type <- 1
+  } else {
+    data$obs_type <- data[[obs_type_label]]
+  }
   colnames(data) <- tolower(colnames(data))
   if(!all(tolower(unlist(cols)) %in% names(data))) {
     stop("Expected column names were not found in data. Please use 'cols' argument to specify column names for independent and dependent variable.")
@@ -18,12 +24,6 @@ parse_input_data <- function(
       data$evid <- 0
     }
   }
-  ## Parse data to include label for obs_type
-  if(is.null(obs_type_label)) {
-    data$obs_type <- 1
-  } else {
-    data$obs_type <- data[[obs_type_label]]
-  }
   if("evid" %in% colnames(data)) {
     data <- data[data$evid == 0,]
   }
diff --git a/R/parse_weights.R b/R/parse_weights.R
index 1c5dbcd..ed76acc 100644
--- a/R/parse_weights.R
+++ b/R/parse_weights.R
@@ -1,6 +1,7 @@
 #' Parse weights argument
 #' 
 #' @inheritParams get_map_estimates
+#' @param t_obs vector of observation times
 #' 
 parse_weights <- function(weights, t_obs) {
   if(!is.null(weights)) {
diff --git a/man/calc_residuals.Rd b/man/calc_residuals.Rd
index 002eb5a..f3c5619 100644
--- a/man/calc_residuals.Rd
+++ b/man/calc_residuals.Rd
@@ -12,6 +12,7 @@ calc_residuals(
   parameters_population,
   covariates,
   regimen,
+  lagtime,
   omega_full,
   error,
   weights,
@@ -42,6 +43,9 @@ calc_residuals(
 
 \item{regimen}{regimen}
 
+\item{lagtime}{vector of lagtimes for each compartment in the model, either
+numeric vector, or vector of parameters.}
+
 \item{omega_full}{full omega matrix}
 
 \item{error}{residual error, specified as list with arguments `add` and/or 
diff --git a/man/get_individual_parameters_from_fit.Rd b/man/get_individual_parameters_from_fit.Rd
index 7fd95bf..9824bbb 100644
--- a/man/get_individual_parameters_from_fit.Rd
+++ b/man/get_individual_parameters_from_fit.Rd
@@ -14,7 +14,11 @@ get_individual_parameters_from_fit(
 \arguments{
 \item{fit}{fit object}
 
-\item{parameters}{list of population parameters}
+\item{parameters}{list of parameters}
+
+\item{nonfixed}{vector of parameters that are not fixed (i.e. estimated)}
+
+\item{as_eta}{vector of parameters that are estimates as eta (e.g. IOV)}
 }
 \description{
 Get the individual parameters from the fitted coefficients (etas)
diff --git a/man/get_map_estimates.Rd b/man/get_map_estimates.Rd
index 02918ab..3a62806 100755
--- a/man/get_map_estimates.Rd
+++ b/man/get_map_estimates.Rd
@@ -23,6 +23,7 @@ get_map_estimates(
   include_omega = TRUE,
   include_error = TRUE,
   regimen = NULL,
+  lagtime = NULL,
   t_init = 0,
   int_step_size = 0.1,
   ll_func = ll_func_PKPDsim,
@@ -93,6 +94,9 @@ estimation.}
 
 \item{regimen}{regimen}
 
+\item{lagtime}{vector of lagtimes for each compartment in the model, either
+numeric vector, or vector of parameters.}
+
 \item{t_init}{initialization time before first dose, default 0.}
 
 \item{int_step_size}{integrator step size passed to PKPDsim}
@@ -102,7 +106,12 @@ in this package.}
 
 \item{optimizer}{optimization library to use, default is `optim`}
 
-\item{method}{optimization method, default `BFGS`}
+\item{method}{optimization method, default `BFGS`. This method is a
+gradient-based method, in which the gradients are computed using finite-
+difference method. When the model contains a step-function where a
+estimated parameter is involved in the step-function (such as in the case of
+lagtime), it is advised to use a non-gradient-based estimation method
+such as `Nelder-Mead` to avoid estimation failures.}
 
 \item{control}{list of options passed to `optim()` function}
 
diff --git a/man/ll_func_PKPDsim.Rd b/man/ll_func_PKPDsim.Rd
index b3aa654..c501cb7 100644
--- a/man/ll_func_PKPDsim.Rd
+++ b/man/ll_func_PKPDsim.Rd
@@ -49,6 +49,7 @@ ll_func_PKPDsim(
   t_init = 0,
   calc_ofv,
   regimen,
+  lagtime = c(0),
   steady_state_analytic,
   include_omega,
   include_error,
@@ -57,13 +58,13 @@ ll_func_PKPDsim(
 )
 }
 \arguments{
-\item{data}{vector of}
+\item{data}{data data.frame with columns `t` and `y` (and possibly evid)}
 
 \item{sim_object}{design (event-table) obtained from PKPDsim to be used in simulations/optimizations}
 
-\item{parameters}{parameter list}
+\item{parameters}{list of parameters}
 
-\item{covariates}{covariates list}
+\item{covariates}{list of covariates, each one created using `PKPDsim::new_coviarate()`}
 
 \item{nonfixed}{non-fixed (i.e. estimated) parameters}
 
@@ -115,11 +116,12 @@ ll_func_PKPDsim(
 
 \item{eta24}{eta24}
 
-\item{error}{error model to use, e.g. `list(add = .5, prop = .15)`}
+\item{error}{residual error, specified as list with arguments `add` and/or 
+`prop` specifying the additive and proportional parts}
 
 \item{transf}{transformation function for data, if needed. E.g. `log(x)`. Default is `function(x) = x`.}
 
-\item{model}{PKPDsim model}
+\item{model}{model, created using `PKPDsim::new_ode_model()`}
 
 \item{omega_full}{full omega matrix}
 
@@ -129,21 +131,28 @@ ll_func_PKPDsim(
 
 \item{sig}{signficance, used in optimization}
 
-\item{weights}{weights for data, generally a vector of weights between 0 and 1.}
+\item{weights}{vector of weights for error. Length of vector should be same 
+as length of observation vector. If NULL (default), all weights are equal. 
+Used in both MAP and NP methods. Note that `weights` argument will also 
+affect residuals (residuals will be scaled too).}
 
-\item{as_eta}{implement as regular eta instead of exponential eta, can be vector.}
+\item{as_eta}{vector of parameters that are estimates as eta (e.g. IOV)}
 
 \item{censoring_idx}{censoring indices, used for <LOQ data.}
 
 \item{censoring_label}{censoring label, used for <LOQ data}
 
-\item{iov_bins}{IOV bins object}
+\item{iov_bins}{bins for inter-occasion variability. Passed unchanged to 
+PKPDsim.}
 
-\item{t_init}{init time for simulations, default 0.}
+\item{t_init}{initialization time before first dose, default 0.}
 
 \item{calc_ofv}{function to calculate OFV based on simulated data and set of parameters and omega matrix}
 
-\item{regimen}{PKPDsim regimen}
+\item{regimen}{regimen}
+
+\item{lagtime}{vector of lagtimes for each compartment in the model, either
+numeric vector, or vector of parameters.}
 
 \item{steady_state_analytic}{list object with settings for steady state MAP estimation.}
 
@@ -151,7 +160,7 @@ ll_func_PKPDsim(
 
 \item{include_error}{include residual error in calculation of OFV?}
 
-\item{verbose}{verbose output?}
+\item{verbose}{show more output}
 
 \item{...}{passed on to PKPDsim}
 }
diff --git a/man/ll_func_generic.Rd b/man/ll_func_generic.Rd
index 5d34f78..103aeb3 100644
--- a/man/ll_func_generic.Rd
+++ b/man/ll_func_generic.Rd
@@ -33,6 +33,7 @@ ll_func_generic(
   parameters,
   covariates = NULL,
   regimen = regimen,
+  lagtime = NULL,
   omega_full = omega_full,
   error = error,
   model,
@@ -42,7 +43,7 @@ ll_func_generic(
 )
 }
 \arguments{
-\item{data}{vector of}
+\item{data}{data data.frame with columns `t` and `y` (and possibly evid)}
 
 \item{eta01}{eta1}
 
@@ -92,21 +93,25 @@ ll_func_generic(
 
 \item{eta24}{eta24}
 
-\item{parameters}{parameter list}
+\item{parameters}{list of parameters}
 
-\item{covariates}{covariates list}
+\item{covariates}{list of covariates, each one created using `PKPDsim::new_coviarate()`}
 
-\item{regimen}{PKPDsim regimen}
+\item{regimen}{regimen}
+
+\item{lagtime}{vector of lagtimes for each compartment in the model, either
+numeric vector, or vector of parameters.}
 
 \item{omega_full}{full omega matrix}
 
-\item{error}{error model to use, e.g. `list(add = .5, prop = .15)`}
+\item{error}{residual error, specified as list with arguments `add` and/or 
+`prop` specifying the additive and proportional parts}
 
-\item{model}{PKPDsim model}
+\item{model}{model, created using `PKPDsim::new_ode_model()`}
 
 \item{sig}{signficance, used in optimization}
 
-\item{verbose}{verbose output?}
+\item{verbose}{show more output}
 
 \item{...}{passed on to PKPDsim}
 }
diff --git a/man/parse_weights.Rd b/man/parse_weights.Rd
index 4450de8..a8ee818 100644
--- a/man/parse_weights.Rd
+++ b/man/parse_weights.Rd
@@ -11,6 +11,8 @@ parse_weights(weights, t_obs)
 as length of observation vector. If NULL (default), all weights are equal. 
 Used in both MAP and NP methods. Note that `weights` argument will also 
 affect residuals (residuals will be scaled too).}
+
+\item{t_obs}{vector of observation times}
 }
 \description{
 Parse weights argument
diff --git a/tests/testthat/setup.R b/tests/testthat/setup.R
new file mode 100644
index 0000000..3fbef06
--- /dev/null
+++ b/tests/testthat/setup.R
@@ -0,0 +1,10 @@
+mod_1cmt_oral_lagtime <- PKPDsim::new_ode_model(
+  code = "
+    dAdt[0] = -KA * A[0]
+    dAdt[1] = +KA * A[0] -(CL/V) * A[1]
+  ",
+  lagtime = c("TLAG", 0),
+  obs = list(cmt = 2, scale = "V"),
+  dose = list(cmt = 1, bioav = 1),
+  parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83)
+)
\ No newline at end of file
diff --git a/tests/testthat/test-get_map_estimates.lagtime.R b/tests/testthat/test-get_map_estimates.lagtime.R
new file mode 100644
index 0000000..3bf5c34
--- /dev/null
+++ b/tests/testthat/test-get_map_estimates.lagtime.R
@@ -0,0 +1,44 @@
+test_that("MAP fit works with oral model with lagtime", {
+  ## Simulate some data
+  reg <- PKPDsim::new_regimen(
+    amt = 100,
+    times = c(0, 24),
+    type = "bolus"
+  )
+  data <- PKPDsim::sim(
+    mod_1cmt_oral_lagtime,
+    regimen = reg,
+    parameters = list(CL = 5, V = 50, KA = 0.5, TLAG = 0.83),
+    t_obs = c(0.25, 0.5, 0.75, 1.0, 2, 5, 8, 16),
+    only_obs = TRUE
+  )
+  ## check that fit works for vector of lagtimes
+  lagtimes <- c(0.25, 0.5, 0.75, 1.0, 2.0)
+  est <- c()
+  for(i in seq(lagtimes)) {
+    fit <- get_map_estimates(
+      model = mod_1cmt_oral_lagtime,
+      parameters = list(CL = 6, V = 40, KA = 0.6, TLAG = lagtimes[i]), # slightly tweaked
+      regimen = reg,
+      omega = c(
+        0.1, 
+        0, 0.1,
+        0, 0, 0.1
+      ),
+      fixed = c("KA"),
+      error = list(prop = 0.1, add = 0.15),
+      data = data,
+      ## !important! due to lagtime discontinuity, cannot use gradient-based optimizers such as
+      ## BGFS or L-BFGS-B, eventhough gradients are based on difference method. Optimization
+      ## with these methods sometimes still work if lagtime initial estimate is higher than
+      ## true value, but better to use a gradient-free method such as Nelder-Mead. 
+      method = "Nelder-Mead",
+      verbose = F
+    )
+    est <- c(est, fit$parameters$TLAG)
+  }
+  expect_equal(
+    round(est, 3),
+    c(0.594, 0.757, 0.856, 1.043, 1.307)
+  )
+})
diff --git a/tests/testthat/test-parse_input_data.R b/tests/testthat/test-parse_input_data.R
index 2457a2f..05b279e 100644
--- a/tests/testthat/test-parse_input_data.R
+++ b/tests/testthat/test-parse_input_data.R
@@ -9,11 +9,11 @@ test_that("parse_input_data handles regular data frame input", {
   result <- parse_input_data(data)
   
   # Check column names are lowercase
-  expect_equal(names(result), c("t", "obs_type", "y"))
+  expect_equal(names(result), c("t", "obs_type", "y", "obs_type"))
   
   # Check sorting
   expect_equal(result$t, c(1, 2, 3))
-  expect_equal(result$obs_type, c(1, 1, 1))
+  expect_equal(result$obs_type, c(1, 2, 1))
 })
 
 test_that("parse_input_data handles PKPDsim object", {
@@ -133,3 +133,4 @@ test_that("parse_input_data preserves data values after transformations", {
   expect_equal(result$dv, c(10, 30))
   expect_equal(result$t, c(2, 3))
 })
+