Merge branch 'main' into custom_lsoda_call

11 files changed, 164 insertions, 37 deletions

diff --git a/man/ds_mixed.Rd b/man/ds_mixed.Rd
new file mode 100644
index 00000000..227b8e7f
--- /dev/null
+++ b/man/ds_mixed.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ds_mixed.R
+\name{ds_mixed}
+\alias{ds_mixed}
+\alias{ds_sfo}
+\alias{ds_fomc}
+\alias{ds_dfop}
+\alias{ds_hs}
+\alias{ds_dfop_sfo}
+\title{Synthetic data for hierarchical kinetic degradation models}
+\description{
+The R code used to create this data object is installed with this package in
+the 'dataset_generation' directory.
+}
+\examples{
+\dontrun{
+  sfo_mmkin <- mmkin("SFO", ds_sfo, quiet = TRUE, error_model = "tc", cores = 15)
+  sfo_saem <- saem(sfo_mmkin, no_random_effect = "parent_0")
+  plot(sfo_saem)
+}
+
+# This is the code used to generate the datasets
+cat(readLines(system.file("dataset_generation/ds_mixed.R", package = "mkin")), sep = "\n")
+}
diff --git a/man/hierarchical_kinetics.Rd b/man/hierarchical_kinetics.Rd
new file mode 100644
index 00000000..2a8e211c
--- /dev/null
+++ b/man/hierarchical_kinetics.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/hierarchical_kinetics.R
+\name{hierarchical_kinetics}
+\alias{hierarchical_kinetics}
+\title{Hierarchical kinetics template}
+\usage{
+hierarchical_kinetics(..., keep_tex = FALSE)
+}
+\arguments{
+\item{...}{Arguments to \code{rmarkdown::pdf_document}}
+
+\item{keep_tex}{Keep the intermediate tex file used in the conversion to PDF}
+}
+\value{
+R Markdown output format to pass to
+\code{\link[rmarkdown:render]{render}}
+}
+\description{
+R markdown format for setting up hierarchical kinetics based on a template
+provided with the mkin package.
+}
+\examples{
+
+\dontrun{
+library(rmarkdown)
+draft("example_analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+}
+
+}
diff --git a/man/illparms.Rd b/man/illparms.Rd
index 14be9c35..75eb18f0 100644
--- a/man/illparms.Rd
+++ b/man/illparms.Rd
@@ -22,7 +22,14 @@ illparms(object, ...)
 
 \method{print}{illparms.mmkin}(x, ...)
 
-\method{illparms}{saem.mmkin}(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...)
+\method{illparms}{saem.mmkin}(
+  object,
+  conf.level = 0.95,
+  random = TRUE,
+  errmod = TRUE,
+  slopes = TRUE,
+  ...
+)
 
 \method{print}{illparms.saem.mmkin}(x, ...)
 
@@ -43,6 +50,10 @@ illparms(object, ...)
 
 \item{errmod}{For hierarchical fits, should error model parameters be
 tested?}
+
+\item{slopes}{For hierarchical \link{saem} fits using saemix as backend,
+should slope parameters in the covariate model(starting with 'beta_') be
+tested?}
 }
 \value{
 For \link{mkinfit} or \link{saem} objects, a character vector of parameter
diff --git a/man/mhmkin.Rd b/man/mhmkin.Rd
index 597fa8ac..c77f4289 100644
--- a/man/mhmkin.Rd
+++ b/man/mhmkin.Rd
@@ -18,7 +18,6 @@ mhmkin(objects, ...)
   backend = "saemix",
   algorithm = "saem",
   no_random_effect = NULL,
-  auto_ranef_threshold = 3,
   ...,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(),
   cluster = NULL
@@ -42,14 +41,14 @@ supported}
 
 \item{algorithm}{The algorithm to be used for fitting (currently not used)}
 
-\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If
-you specify "auto", random effects are only included if the number
-of datasets in which the parameter passed the t-test is at least 'auto_ranef_threshold'.
-Beware that while this may make for convenient model reduction or even
-numerical stability of the algorithm, it will likely lead to
-underparameterised models.}
-
-\item{auto_ranef_threshold}{See 'no_random_effect.}
+\item{no_random_effect}{Default is NULL and will be passed to \link{saem}. If a
+character vector is supplied, it will be passed to all calls to \link{saem},
+which will exclude random effects for all matching parameters. Alternatively,
+a list of character vectors or an object of class \link{illparms.mhmkin} can be
+specified. They have to have the same dimensions that the return object of
+the current call will have, i.e. the number of rows must match the number
+of degradation models in the mmkin object(s), and the number of columns must
+match the number of error models used in the mmkin object(s).}
 
 \item{cores}{The number of cores to be used for multicore processing. This
 is only used when the \code{cluster} argument is \code{NULL}. On Windows
@@ -76,7 +75,7 @@ be indexed using the degradation model names for the first index (row index)
 and the error model names for the second index (column index), with class
 attribute 'mhmkin'.
 
-An object of class \code{\link{mhmkin}}.
+An object inheriting from \code{\link{mhmkin}}.
 }
 \description{
 The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchichal
@@ -84,6 +83,43 @@ The name of the methods expresses that (\strong{m}ultiple) \strong{h}ierarchicha
 fitted. Our kinetic models are nonlinear, so we can use various nonlinear
 mixed-effects model fitting functions.
 }
+\examples{
+\dontrun{
+# We start with separate evaluations of all the first six datasets with two
+# degradation models and two error models
+f_sep_const <- mmkin(c("SFO", "FOMC"), ds_fomc[1:6], cores = 2, quiet = TRUE)
+f_sep_tc <- update(f_sep_const, error_model = "tc")
+# The mhmkin function sets up hierarchical degradation models aka
+# nonlinear mixed-effects models for all four combinations, specifying
+# uncorrelated random effects for all degradation parameters
+f_saem_1 <- mhmkin(list(f_sep_const, f_sep_tc), cores = 2)
+status(f_saem_1)
+# The 'illparms' function shows that in all hierarchical fits, at least
+# one random effect is ill-defined (the confidence interval for the
+# random effect expressed as standard deviation includes zero)
+illparms(f_saem_1)
+# Therefore we repeat the fits, excluding the ill-defined random effects
+f_saem_2 <- update(f_saem_1, no_random_effect = illparms(f_saem_1))
+status(f_saem_2)
+illparms(f_saem_2)
+# Model comparisons show that FOMC with two-component error is preferable,
+# and confirms our reduction of the default parameter model
+anova(f_saem_1)
+anova(f_saem_2)
+# The convergence plot for the selected model looks fine
+saemix::plot(f_saem_2[["FOMC", "tc"]]$so, plot.type = "convergence")
+# The plot of predictions versus data shows that we have a pretty data-rich
+# situation with homogeneous distribution of residuals, because we used the
+# same degradation model, error model and parameter distribution model that
+# was used in the data generation.
+plot(f_saem_2[["FOMC", "tc"]])
+# We can specify the same parameter model reductions manually
+no_ranef <- list("parent_0", "log_beta", "parent_0", c("parent_0", "log_beta"))
+dim(no_ranef) <- c(2, 2)
+f_saem_2m <- update(f_saem_1, no_random_effect = no_ranef)
+anova(f_saem_2m)
+}
+}
 \seealso{
 \code{\link{[.mhmkin}} for subsetting \link{mhmkin} objects
 }
diff --git a/man/mkinmod.Rd b/man/mkinmod.Rd
index 612c3c2b..77a4f520 100644
--- a/man/mkinmod.Rd
+++ b/man/mkinmod.Rd
@@ -33,7 +33,7 @@ the source compartment.
 Additionally, \code{\link[=mkinsub]{mkinsub()}} has an argument \code{to}, specifying names of
 variables to which a transfer is to be assumed in the model.
 If the argument \code{use_of_ff} is set to "min"
-(default) and the model for the compartment is "SFO" or "SFORB", an
+and the model for the compartment is "SFO" or "SFORB", an
 additional \code{\link[=mkinsub]{mkinsub()}} argument can be \code{sink = FALSE}, effectively
 fixing the flux to sink to zero.
 In print.mkinmod, this argument is currently not used.}
diff --git a/man/parplot.Rd b/man/parplot.Rd
index 37c5841d..67bf0cc1 100644
--- a/man/parplot.Rd
+++ b/man/parplot.Rd
@@ -10,6 +10,7 @@ parplot(object, ...)
 \method{parplot}{multistart.saem.mmkin}(
   object,
   llmin = -Inf,
+  llquant = NA,
   scale = c("best", "median"),
   lpos = "bottomleft",
   main = "",
@@ -23,7 +24,11 @@ parplot(object, ...)
 
 \item{llmin}{The minimum likelihood of objects to be shown}
 
-\item{scale}{By default, scale parameters using the best available fit.
+\item{llquant}{Fractional value for selecting only the fits with higher
+likelihoods. Overrides 'llmin'.}
+
+\item{scale}{By default, scale parameters using the best
+available fit.
 If 'median', parameters are scaled using the median parameters from all fits.}
 
 \item{lpos}{Positioning of the legend.}
@@ -35,6 +40,11 @@ Produces a boxplot with all parameters from the multiple runs, scaled
 either by the parameters of the run with the highest likelihood,
 or by their medians as proposed in the paper by Duchesne et al. (2021).
 }
+\details{
+Starting values of degradation model parameters and error model parameters
+are shown as green circles. The results obtained in the original run
+are shown as red circles.
+}
 \references{
 Duchesne R, Guillemin A, Gandrillon O, Crauste F. Practical
 identifiability in the frame of nonlinear mixed effects models: the example
diff --git a/man/read_spreadsheet.Rd b/man/read_spreadsheet.Rd
index 147d09bf..41c32108 100644
--- a/man/read_spreadsheet.Rd
+++ b/man/read_spreadsheet.Rd
@@ -7,7 +7,7 @@
 read_spreadsheet(
   path,
   valid_datasets = "all",
-  parent_only = TRUE,
+  parent_only = FALSE,
   normalize = TRUE
 )
 }
diff --git a/man/saem.Rd b/man/saem.Rd
index 3a5abada..89647ff4 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -24,7 +24,7 @@ saem(object, ...)
   covariates = NULL,
   covariate_models = NULL,
   no_random_effect = NULL,
-  error.init = c(3, 0.1),
+  error.init = c(1, 1),
   nbiter.saemix = c(300, 100),
   control = list(displayProgress = FALSE, print = FALSE, nbiter.saemix = nbiter.saemix,
     save = FALSE, save.graphs = FALSE),
diff --git a/man/summary.saem.mmkin.Rd b/man/summary.saem.mmkin.Rd
index fb099899..0845d4d2 100644
--- a/man/summary.saem.mmkin.Rd
+++ b/man/summary.saem.mmkin.Rd
@@ -92,10 +92,21 @@ f_mmkin_dfop_sfo <- mmkin(list(dfop_sfo), ds_syn_dfop_sfo,
 f_saem_dfop_sfo <- saem(f_mmkin_dfop_sfo)
 print(f_saem_dfop_sfo)
 illparms(f_saem_dfop_sfo)
-f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo, covariance.model = diag(c(0, 0, 1, 1, 1, 0)))
+f_saem_dfop_sfo_2 <- update(f_saem_dfop_sfo,
+  no_random_effect = c("parent_0", "log_k_m1"))
 illparms(f_saem_dfop_sfo_2)
 intervals(f_saem_dfop_sfo_2)
 summary(f_saem_dfop_sfo_2, data = TRUE)
+# Add a correlation between random effects of g and k2
+cov_model_3 <- f_saem_dfop_sfo_2$so@model@covariance.model
+cov_model_3["log_k2", "g_qlogis"] <- 1
+cov_model_3["g_qlogis", "log_k2"] <- 1
+f_saem_dfop_sfo_3 <- update(f_saem_dfop_sfo,
+  covariance.model = cov_model_3)
+intervals(f_saem_dfop_sfo_3)
+# The correlation does not improve the fit judged by AIC and BIC, although
+# the likelihood is higher with the additional parameter
+anova(f_saem_dfop_sfo, f_saem_dfop_sfo_2, f_saem_dfop_sfo_3)
 }
 
 }
diff --git a/man/summary_listing.Rd b/man/summary_listing.Rd
new file mode 100644
index 00000000..995ebd8d
--- /dev/null
+++ b/man/summary_listing.Rd
@@ -0,0 +1,27 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/summary_listing.R
+\name{summary_listing}
+\alias{summary_listing}
+\alias{tex_listing}
+\alias{html_listing}
+\title{Display the output of a summary function according to the output format}
+\usage{
+summary_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+tex_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
+
+html_listing(object, caption = NULL)
+}
+\arguments{
+\item{object}{The object for which the summary is to be listed}
+
+\item{caption}{An optional caption}
+
+\item{label}{An optional label, ignored in html output}
+
+\item{clearpage}{Should a new page be started after the listing? Ignored in html output}
+}
+\description{
+This function is intended for use in a R markdown code chunk with the chunk
+option \code{results = "asis"}.
+}
diff --git a/man/tex_listing.Rd b/man/tex_listing.Rd
deleted file mode 100644
index 2f11d211..00000000
--- a/man/tex_listing.Rd
+++ /dev/null
@@ -1,21 +0,0 @@
-% Generated by roxygen2: do not edit by hand
-% Please edit documentation in R/tex_listing.R
-\name{tex_listing}
-\alias{tex_listing}
-\title{Wrap the output of a summary function in tex listing environment}
-\usage{
-tex_listing(object, caption = NULL, label = NULL, clearpage = TRUE)
-}
-\arguments{
-\item{object}{The object for which the summary is to be listed}
-
-\item{caption}{An optional caption}
-
-\item{label}{An optional label}
-
-\item{clearpage}{Should a new page be started after the listing?}
-}
-\description{
-This function can be used in a R markdown code chunk with the chunk
-option \code{results = "asis"}.
-}