Complete documentation and fix a bug

The bug was introduced by the changes in summary.saem.mmkin.R and
surfaced in the tests when using saemix transformations.

11 files changed, 69 insertions, 29 deletions

diff --git a/NAMESPACE b/NAMESPACE
index eb34304e..aaaa6484 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -147,6 +147,7 @@ importFrom(KernSmooth,bkde)
 importFrom(R6,R6Class)
 importFrom(grDevices,dev.cur)
 importFrom(lmtest,lrtest)
+importFrom(methods,is)
 importFrom(methods,signature)
 importFrom(nlme,intervals)
 importFrom(parallel,detectCores)
@@ -170,6 +171,7 @@ importFrom(stats,na.fail)
 importFrom(stats,nlminb)
 importFrom(stats,nobs)
 importFrom(stats,optimize)
+importFrom(stats,pchisq)
 importFrom(stats,plogis)
 importFrom(stats,predict)
 importFrom(stats,pt)
@@ -181,8 +183,10 @@ importFrom(stats,qt)
 importFrom(stats,residuals)
 importFrom(stats,rnorm)
 importFrom(stats,shapiro.test)
+importFrom(stats,terms)
 importFrom(stats,update)
 importFrom(stats,vcov)
+importFrom(utils,capture.output)
 importFrom(utils,getFromNamespace)
 importFrom(utils,packageVersion)
 importFrom(utils,write.table)
diff --git a/R/anova.saem.mmkin.R b/R/anova.saem.mmkin.R
index e88561b9..9937a919 100644
--- a/R/anova.saem.mmkin.R
+++ b/R/anova.saem.mmkin.R
@@ -13,7 +13,9 @@
 #' the alternative models are tested against the first model. Should
 #' only be done for nested models.
 #' @param model.names Optional character vector of model names
-#' @importFrom stats anova logLik update
+#' @importFrom stats anova logLik update pchisq terms
+#' @importFrom methods is
+#' @importFrom utils capture.output
 #' @export
 #' @return an "anova" data frame; the traditional (S3) result of anova()
 anova.saem.mmkin <- function(object, ...,
diff --git a/R/saem.R b/R/saem.R
index 4b94a3d2..79b4b9ee 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -723,6 +723,7 @@ saemix_data <- function(object, covariates = NULL, verbose = FALSE, ...) {
 #' logLik method for saem.mmkin objects
 #'
 #' @param object The fitted [saem.mmkin] object
+#' @param \dots Passed to [saemix::logLik.SaemixObject]
 #' @param method Passed to [saemix::logLik.SaemixObject]
 #' @export
 logLik.saem.mmkin <- function(object, ..., method = c("lin", "is", "gq")) {
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index c66294ab..651cc29e 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -132,7 +132,7 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
       }
     }
   } else {
-    confint_back <- confint_trans
+    confint_back <- confint_trans[names_fixed_effects, ]
   }
 
   #  Correlation of fixed effects (inspired by summary.nlme)
@@ -210,6 +210,7 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
   cat("Using", paste(x$so@options$nbiter.saemix, collapse = ", "),
     "iterations and", x$so@options$nb.chains, "chains\n")
 
+  cat("\nVariance model: ")
   cat(switch(x$err_mod,
     const = "Constant variance",
     obs = "Variance unique to each observed variable",
diff --git a/log/check.log b/log/check.log
index 8fce3fd1..c2b4aa50 100644
--- a/log/check.log
+++ b/log/check.log
@@ -7,8 +7,10 @@
 * checking extension type ... Package
 * this is package ‘mkin’ version ‘1.1.2’
 * package encoding: UTF-8
-* checking CRAN incoming feasibility ... Note_to_CRAN_maintainers
+* checking CRAN incoming feasibility ... NOTE
 Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
+
+The Date field is over a month old.
 * checking package namespace information ... OK
 * checking package dependencies ... OK
 * checking if this is a source package ... OK
@@ -69,5 +71,9 @@ Maintainer: ‘Johannes Ranke <johannes.ranke@jrwb.de>’
 * checking for detritus in the temp directory ... OK
 * DONE
 
-Status: OK
+Status: 1 NOTE
+See
+  ‘/home/jranke/git/mkin/mkin.Rcheck/00check.log’
+for details.
+
 
diff --git a/log/test.log b/log/test.log
index d7d623d7..b6ec2084 100644
--- a/log/test.log
+++ b/log/test.log
@@ -1,11 +1,11 @@
 ℹ Testing mkin
 ✔ | F W S  OK | Context
 ✔ |         5 | AIC calculation
-✔ |         5 | Analytical solutions for coupled models [3.2s]
+✔ |         5 | Analytical solutions for coupled models [3.3s]
 ✔ |         5 | Calculation of Akaike weights
 ✔ |         3 | Export dataset for reading into CAKE
 ✔ |        12 | Confidence intervals and p-values [1.0s]
-✔ |     1  12 | Dimethenamid data from 2018 [31.3s]
+✔ |     1  12 | Dimethenamid data from 2018 [31.6s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_dmta.R:98:3): Different backends get consistent results for SFO-SFO3+, dimethenamid data
 Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
@@ -16,7 +16,7 @@ Reason: Fitting this ODE model with saemix takes about 15 minutes on my system
 ✔ |        14 | Results for FOCUS D established in expertise for UBA (Ranke 2014) [0.8s]
 ✔ |         4 | Test fitting the decline of metabolites from their maximum [0.4s]
 ✔ |         1 | Fitting the logistic model [0.2s]
-✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [14.4s]
+✔ |         7 | Batch fitting and diagnosing hierarchical kinetic models [14.3s]
 ✔ |     1  12 | Nonlinear mixed-effects models [0.3s]
 ────────────────────────────────────────────────────────────────────────────────
 Skip (test_mixed.R:68:3): saemix results are reproducible for biphasic fits
@@ -29,20 +29,20 @@ Reason: Fitting with saemix takes around 10 minutes when using deSolve
 ✔ |         3 | Model predictions with mkinpredict [0.4s]
 ✔ |        16 | Evaluations according to 2015 NAFTA guidance [1.8s]
 ✔ |         9 | Nonlinear mixed-effects models with nlme [8.3s]
-✔ |        16 | Plotting [9.7s]
+✔ |        16 | Plotting [9.9s]
 ✔ |         4 | Residuals extracted from mkinfit models
-✔ |        37 | saemix parent models [208.2s]
-✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.4s]
+✔ |        37 | saemix parent models [208.7s]
+✔ |         2 | Complex test case from Schaefer et al. (2007) Piacenza paper [1.3s]
 ✔ |        11 | Processing of residue series
-✔ |         7 | Fitting the SFORB model [3.6s]
+✔ |         7 | Fitting the SFORB model [3.8s]
 ✔ |         1 | Summaries of old mkinfit objects
 ✔ |         5 | Summary [0.2s]
-✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.1s]
-✔ |         9 | Hypothesis tests [7.6s]
-✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.1s]
+✔ |         4 | Results for synthetic data established in expertise for UBA (Ranke 2014) [2.2s]
+✔ |         9 | Hypothesis tests [7.8s]
+✔ |         4 | Calculation of maximum time weighted average concentrations (TWAs) [2.2s]
 
 ══ Results ═════════════════════════════════════════════════════════════════════
-Duration: 304.2 s
+Duration: 305.6 s
 
 ── Skipped tests  ──────────────────────────────────────────────────────────────
 • Fitting this ODE model with saemix takes about 15 minutes on my system (1)
diff --git a/man/anova.saem.mmkin.Rd b/man/anova.saem.mmkin.Rd
index 8f9cb241..ab6022bc 100644
--- a/man/anova.saem.mmkin.Rd
+++ b/man/anova.saem.mmkin.Rd
@@ -31,5 +31,6 @@ only be done for nested models.}
 an "anova" data frame; the traditional (S3) result of anova()
 }
 \description{
-Anova method for saem.mmkin objects
+Generate an anova object. The method to calculate the BIC is that from
+the saemix package. As in other prominent anova methods, models are sorted
 }
diff --git a/man/logLik.saem.mmkin.Rd b/man/logLik.saem.mmkin.Rd
index 8ea3426f..603f4607 100644
--- a/man/logLik.saem.mmkin.Rd
+++ b/man/logLik.saem.mmkin.Rd
@@ -7,6 +7,10 @@
 \method{logLik}{saem.mmkin}(object, ..., method = c("lin", "is", "gq"))
 }
 \arguments{
+\item{object}{The fitted \link{saem.mmkin} object}
+
+\item{\dots}{Passed to \link[saemix:logLik]{saemix::logLik.SaemixObject}}
+
 \item{method}{Passed to \link[saemix:logLik]{saemix::logLik.SaemixObject}}
 }
 \description{
diff --git a/man/saem.Rd b/man/saem.Rd
index 2b9199dd..d7b04691 100644
--- a/man/saem.Rd
+++ b/man/saem.Rd
@@ -82,6 +82,13 @@ automatic choice is not desired}
 default, uncorrelated random effects are specified for all degradation
 parameters.}
 
+\item{covariates}{A data frame with covariate data for use in
+'covariate_models', with dataset names as row names.}
+
+\item{covariate_models}{A list containing linear model formulas with one explanatory
+variable, i.e. of the type 'parameter ~ covariate'. Covariates must be available
+in the 'covariates' data frame.}
+
 \item{no_random_effect}{Character vector of degradation parameters for
 which there should be no variability over the groups. Only used
 if the covariance model is not explicitly specified.}
@@ -142,10 +149,13 @@ f_mmkin_parent <- mmkin(c("SFO", "FOMC", "DFOP"), ds, quiet = TRUE)
 f_saem_sfo <- saem(f_mmkin_parent["SFO", ])
 f_saem_fomc <- saem(f_mmkin_parent["FOMC", ])
 f_saem_dfop <- saem(f_mmkin_parent["DFOP", ])
+anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
+anova(f_saem_sfo, f_saem_dfop, test = TRUE)
 illparms(f_saem_dfop)
-update(f_saem_dfop, covariance.model = diag(c(1, 1, 1, 0)))
-AIC(f_saem_dfop)
+f_saem_dfop_red <- update(f_saem_dfop, no_random_effect = "g_qlogis")
+anova(f_saem_dfop, f_saem_dfop_red, test = TRUE)
 
+anova(f_saem_sfo, f_saem_fomc, f_saem_dfop)
 # The returned saem.mmkin object contains an SaemixObject, therefore we can use
 # functions from saemix
 library(saemix)
@@ -157,7 +167,7 @@ plot(f_saem_fomc$so, plot.type = "vpc")
 
 f_mmkin_parent_tc <- update(f_mmkin_parent, error_model = "tc")
 f_saem_fomc_tc <- saem(f_mmkin_parent_tc["FOMC", ])
-compare.saemix(f_saem_fomc$so, f_saem_fomc_tc$so)
+anova(f_saem_fomc, f_saem_fomc_tc, test = TRUE)
 
 sfo_sfo <- mkinmod(parent = mkinsub("SFO", "A1"),
   A1 = mkinsub("SFO"))
diff --git a/tests/testthat/summary_hfit_sfo_tc.txt b/tests/testthat/summary_hfit_sfo_tc.txt
index e3e2f7e4..49765187 100644
--- a/tests/testthat/summary_hfit_sfo_tc.txt
+++ b/tests/testthat/summary_hfit_sfo_tc.txt
@@ -31,9 +31,12 @@ Likelihood computed by importance sampling
   533 531   -261
 
 Optimised parameters:
-             est. lower upper
-parent_0      101   100   102
-log_k_parent   -3    -4    -3
+                  est.  lower  upper
+parent_0        101.01  99.57 102.45
+log_k_parent     -3.32  -3.53  -3.11
+a.1               0.90   0.64   1.17
+b.1               0.05   0.04   0.06
+SD.log_k_parent   0.27   0.11   0.42
 
 Correlation: 
              pr_0
diff --git a/tests/testthat/summary_saem_biphasic_s.txt b/tests/testthat/summary_saem_biphasic_s.txt
index 6b203991..7c337843 100644
--- a/tests/testthat/summary_saem_biphasic_s.txt
+++ b/tests/testthat/summary_saem_biphasic_s.txt
@@ -38,13 +38,21 @@ Likelihood computed by importance sampling
   2334 2344  -1153
 
 Optimised parameters:
-                est. lower upper
-parent_0       1e+02 1e+02 1e+02
-k_m1           7e-03 6e-03 7e-03
-f_parent_to_m1 5e-01 4e-01 5e-01
-k1             1e-01 9e-02 1e-01
-k2             2e-02 2e-02 3e-02
-g              5e-01 5e-01 5e-01
+                   est.  lower upper
+parent_0          1e+02  1e+02 1e+02
+k_m1              7e-03  6e-03 7e-03
+f_parent_to_m1    5e-01  4e-01 5e-01
+k1                1e-01  9e-02 1e-01
+k2                2e-02  2e-02 3e-02
+g                 5e-01  5e-01 5e-01
+a.1               9e-01  8e-01 1e+00
+b.1               5e-02  5e-02 6e-02
+SD.parent_0       3e-02 -5e+01 5e+01
+SD.k_m1           2e-01  1e-01 3e-01
+SD.f_parent_to_m1 3e-01  2e-01 4e-01
+SD.k1             4e-01  2e-01 5e-01
+SD.k2             3e-01  2e-01 5e-01
+SD.g              2e-01  6e-02 4e-01
 
 Correlation: 
                pr_0 k_m1 f___ k1   k2