diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/confint.mkinfit.R | 8 | ||||
| -rw-r--r-- | R/logLik.mkinfit.R | 6 | ||||
| -rw-r--r-- | R/mkinfit.R | 2 | ||||
| -rw-r--r-- | R/mkinresplot.R | 2 | ||||
| -rw-r--r-- | R/mmkin.R | 1 | ||||
| -rw-r--r-- | R/nlme.mmkin.R | 36 | ||||
| -rw-r--r-- | R/transform_odeparms.R | 18 |
7 files changed, 40 insertions, 33 deletions
diff --git a/R/confint.mkinfit.R b/R/confint.mkinfit.R index 53eb45ee..6403c349 100644 --- a/R/confint.mkinfit.R +++ b/R/confint.mkinfit.R @@ -57,20 +57,20 @@ #' if (identical(Sys.getenv("NOT_CRAN"), "true")) { #' n_cores <- parallel::detectCores() - 1 #' } else { -#' n_cores <- 1 +#' n_cores <- 1 #' } #' if (Sys.getenv("TRAVIS") != "") n_cores = 1 #' if (Sys.info()["sysname"] == "Windows") n_cores = 1 #' -#' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), quiet = TRUE) +#' SFO_SFO <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), +#' use_of_ff = "min", quiet = TRUE) #' SFO_SFO.ff <- mkinmod(parent = mkinsub("SFO", "m1"), m1 = mkinsub("SFO"), #' use_of_ff = "max", quiet = TRUE) #' f_d_1 <- mkinfit(SFO_SFO, subset(FOCUS_2006_D, value != 0), quiet = TRUE) #' system.time(ci_profile <- confint(f_d_1, method = "profile", cores = 1, quiet = TRUE)) #' # Using more cores does not save much time here, as parent_0 takes up most of the time #' # If we additionally exclude parent_0 (the confidence of which is often of -#' # minor interest), we get a nice performance improvement from about 50 -#' # seconds to about 12 seconds if we use at least four cores +#' # minor interest), we get a nice performance improvement if we use at least 4 cores #' system.time(ci_profile_no_parent_0 <- confint(f_d_1, method = "profile", #' c("k_parent_sink", "k_parent_m1", "k_m1_sink", "sigma"), cores = n_cores)) #' ci_profile diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R index 1c025893..7cc10234 100644 --- a/R/logLik.mkinfit.R +++ b/R/logLik.mkinfit.R @@ -25,10 +25,10 @@ #' parent = mkinsub("SFO", to = "m1"), #' m1 = mkinsub("SFO") #' ) -#' d_t <- FOCUS_2006_D +#' d_t <- subset(FOCUS_2006_D, value != 0) #' f_nw <- mkinfit(sfo_sfo, d_t, quiet = TRUE) # no weighting (weights are unity) -#' f_obs <- mkinfit(sfo_sfo, d_t, error_model = "obs", quiet = TRUE) -#' f_tc <- mkinfit(sfo_sfo, d_t, error_model = "tc", quiet = TRUE) +#' f_obs <- update(f_nw, error_model = "obs") +#' f_tc <- update(f_nw, error_model = "tc") #' AIC(f_nw, f_obs, f_tc) #' } #' diff --git a/R/mkinfit.R b/R/mkinfit.R index d7de718b..704e70a9 100644 --- a/R/mkinfit.R +++ b/R/mkinfit.R @@ -542,6 +542,8 @@ mkinfit <- function(mkinmod, observed, { assign("calls", calls + 1, inherits = TRUE) # Increase the model solution counter + #browser() + # Trace parameter values if requested and if we are actually optimising if(trace_parms & update_data) cat(format(P, width = 10, digits = 6), "\n") diff --git a/R/mkinresplot.R b/R/mkinresplot.R index bad28ae8..be361690 100644 --- a/R/mkinresplot.R +++ b/R/mkinresplot.R @@ -28,7 +28,7 @@ utils::globalVariables(c("variable", "residual")) #' @param \dots further arguments passed to \code{\link{plot}}. #' @return Nothing is returned by this function, as it is called for its side #' effect, namely to produce a plot. -#' @author Johannes Ranke +#' @author Johannes Ranke and Katrin Lindenberger #' @seealso \code{\link{mkinplot}}, for a way to plot the data and the fitted #' lines of the mkinfit object, and \code{\link{plot_res}} for a function #' combining the plot of the fit and the residual plot. @@ -162,6 +162,7 @@ mmkin <- function(models = c("SFO", "FOMC", "DFOP"), datasets, #' #' @param x An [mmkin] object. #' @param \dots Not used. +#' @rdname mmkin #' @export print.mmkin <- function(x, ...) { cat("<mmkin> object\n") diff --git a/R/nlme.mmkin.R b/R/nlme.mmkin.R index 82d5f6de..ff1f2fff 100644 --- a/R/nlme.mmkin.R +++ b/R/nlme.mmkin.R @@ -24,6 +24,11 @@ get_deg_func <- function() { #' This functions sets up a nonlinear mixed effects model for an mmkin row #' object. An mmkin row object is essentially a list of mkinfit objects that #' have been obtained by fitting the same model to a list of datasets. +#' +#' Note that the convergence of the nlme algorithms depends on the quality +#' of the data. In degradation kinetics, we often only have few datasets +#' (e.g. data for few soils) and complicated degradation models, which may +#' make it impossible to obtain convergence with nlme. #' #' @param model An [mmkin] row object. #' @param data Ignored, data are taken from the mmkin model @@ -88,11 +93,10 @@ get_deg_func <- function() { #' # f_nlme_sfo_sfo_ff <- nlme(f_2["SFO-SFO-ff", ]) #' #plot(f_nlme_sfo_sfo_ff) #' -#' # With the log-Cholesky parameterization, this converges in 11 -#' # iterations and around 100 seconds, but without tweaking control -#' # parameters (with pdDiag, increasing the tolerance and pnlsMaxIter was -#' # necessary) -#' f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ]) +#' # For the following, we need to increase pnlsMaxIter and the tolerance +#' # to get convergence +#' f_nlme_dfop_sfo <- nlme(f_2["DFOP-SFO", ], +#' control = list(pnlsMaxIter = 120, tolerance = 5e-4)) #' #' plot(f_nlme_dfop_sfo) #' @@ -112,22 +116,18 @@ get_deg_func <- function() { #' print(f_nlme_dfop_tc) #' } #' -#' f_2_obs <- mmkin(list("SFO-SFO" = m_sfo_sfo, -#' "DFOP-SFO" = m_dfop_sfo), -#' ds_2, quiet = TRUE, error_model = "obs") +#' f_2_obs <- update(f_2, error_model = "obs") #' f_nlme_sfo_sfo_obs <- nlme(f_2_obs["SFO-SFO", ]) #' print(f_nlme_sfo_sfo_obs) -#' f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ]) +#' f_nlme_dfop_sfo_obs <- nlme(f_2_obs["DFOP-SFO", ], +#' control = list(pnlsMaxIter = 120, tolerance = 5e-4)) #' -#' f_2_tc <- mmkin(list("SFO-SFO" = m_sfo_sfo, -#' "DFOP-SFO" = m_dfop_sfo), -#' ds_2, quiet = TRUE, error_model = "tc") -#' # f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ]) # stops with error message -#' f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ]) -#' # We get warnings about false convergence in the LME step in several iterations -#' # but as the last such warning occurs in iteration 25 and we have 28 iterations -#' # we can ignore these -#' anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs, f_nlme_dfop_sfo_tc) +#' f_2_tc <- update(f_2, error_model = "tc") +#' # f_nlme_sfo_sfo_tc <- nlme(f_2_tc["SFO-SFO", ]) # No convergence with 50 iterations +#' # f_nlme_dfop_sfo_tc <- nlme(f_2_tc["DFOP-SFO", ], +#' # control = list(pnlsMaxIter = 120, tolerance = 5e-4)) # Error in X[, fmap[[nm]]] <- gradnm +#' +#' anova(f_nlme_dfop_sfo, f_nlme_dfop_sfo_obs) #' #' } nlme.mmkin <- function(model, data = "auto", diff --git a/R/transform_odeparms.R b/R/transform_odeparms.R index f21d31fc..4fe4e5c2 100644 --- a/R/transform_odeparms.R +++ b/R/transform_odeparms.R @@ -42,17 +42,21 @@ #' #' SFO_SFO <- mkinmod( #' parent = list(type = "SFO", to = "m1", sink = TRUE), -#' m1 = list(type = "SFO")) +#' m1 = list(type = "SFO"), use_of_ff = "min") +#' #' # Fit the model to the FOCUS example dataset D using defaults -#' fit <- mkinfit(SFO_SFO, FOCUS_2006_D, quiet = TRUE) +#' FOCUS_D <- subset(FOCUS_2006_D, value != 0) # remove zero values to avoid warning +#' fit <- mkinfit(SFO_SFO, FOCUS_D, quiet = TRUE) #' fit.s <- summary(fit) #' # Transformed and backtransformed parameters #' print(fit.s$par, 3) #' print(fit.s$bpar, 3) #' #' \dontrun{ -#' # Compare to the version without transforming rate parameters -#' fit.2 <- mkinfit(SFO_SFO, FOCUS_2006_D, transform_rates = FALSE, quiet = TRUE) +#' # Compare to the version without transforming rate parameters (does not work +#' # with analytical solution, we get NA values for m1 in predictions) +#' fit.2 <- mkinfit(SFO_SFO, FOCUS_D, transform_rates = FALSE, +#' solution_type = "deSolve", quiet = TRUE) #' fit.2.s <- summary(fit.2) #' print(fit.2.s$par, 3) #' print(fit.2.s$bpar, 3) @@ -66,13 +70,13 @@ #' backtransform_odeparms(transformed, SFO_SFO) #' #' \dontrun{ -#' # The case of formation fractions +#' # The case of formation fractions (this is now the default) #' SFO_SFO.ff <- mkinmod( #' parent = list(type = "SFO", to = "m1", sink = TRUE), #' m1 = list(type = "SFO"), #' use_of_ff = "max") #' -#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_2006_D, quiet = TRUE) +#' fit.ff <- mkinfit(SFO_SFO.ff, FOCUS_D, quiet = TRUE) #' fit.ff.s <- summary(fit.ff) #' print(fit.ff.s$par, 3) #' print(fit.ff.s$bpar, 3) @@ -87,7 +91,7 @@ #' use_of_ff = "max") #' #' -#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_2006_D, quiet = TRUE) +#' fit.ff.2 <- mkinfit(SFO_SFO.ff.2, FOCUS_D, quiet = TRUE) #' fit.ff.2.s <- summary(fit.ff.2) #' print(fit.ff.2.s$par, 3) #' print(fit.ff.2.s$bpar, 3) |