Merge branch 'main' into custom_lsoda_call

15 files changed, 277 insertions, 78 deletions

diff --git a/R/ds_mixed.R b/R/ds_mixed.R
new file mode 100644
index 00000000..c5055712
--- /dev/null
+++ b/R/ds_mixed.R
@@ -0,0 +1,17 @@
+#'  Synthetic data for hierarchical kinetic degradation models
+#'
+#' The R code used to create this data object is installed with this package in
+#' the 'dataset_generation' directory.
+#'
+#' @name ds_mixed
+#' @aliases ds_sfo ds_fomc ds_dfop ds_hs ds_dfop_sfo
+#' @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")
+NULL
diff --git a/R/hierarchical_kinetics.R b/R/hierarchical_kinetics.R
new file mode 100644
index 00000000..e545754a
--- /dev/null
+++ b/R/hierarchical_kinetics.R
@@ -0,0 +1,40 @@
+#' Hierarchical kinetics template
+#'
+#' R markdown format for setting up hierarchical kinetics based on a template
+#' provided with the mkin package.
+#'
+#' @inheritParams rmarkdown::pdf_document
+#' @param ... Arguments to \code{rmarkdown::pdf_document}
+#'
+#' @return R Markdown output format to pass to
+#'   \code{\link[rmarkdown:render]{render}}
+#'
+#' @examples
+#'
+#' \dontrun{
+#' library(rmarkdown)
+#' draft("example_analysis.rmd", template = "hierarchical_kinetics", package = "mkin")
+#' }
+#'
+#' @export
+hierarchical_kinetics <- function(..., keep_tex = FALSE) {
+
+  if (getRversion() < "4.1.0")
+    stop("You need R with version > 4.1.0 to compile this document")
+
+  if (!requireNamespace("knitr")) stop("Please install the knitr package to use this template")
+  if (!requireNamespace("rmarkdown")) stop("Please install the rmarkdown package to use this template")
+  knitr::opts_chunk$set(echo = FALSE, cache = TRUE, comment = "", tidy = FALSE, echo = TRUE)
+  knitr::opts_chunk$set(fig.align = "center", fig.pos = "H")
+  options(knitr.kable.NA = "")
+
+  fmt <- rmarkdown::pdf_document(...,
+    keep_tex = keep_tex,
+    toc = TRUE,
+    toc_depth = 3,
+    includes = rmarkdown::includes(in_header = "header.tex"),
+    extra_dependencies = c("float", "listing", "framed")
+  )
+
+  return(fmt)
+}
diff --git a/R/illparms.R b/R/illparms.R
index 01e75cf1..eef4bd33 100644
--- a/R/illparms.R
+++ b/R/illparms.R
@@ -20,6 +20,9 @@
 #' @param random For hierarchical fits, should random effects be tested?
 #' @param errmod For hierarchical fits, should error model parameters be
 #' tested?
+#' @param slopes For hierarchical [saem] fits using saemix as backend,
+#' should slope parameters in the covariate model(starting with 'beta_') be
+#' tested?
 #' @return For [mkinfit] or [saem] objects, a character vector of parameter
 #' names. For [mmkin] or [mhmkin] objects, a matrix like object of class
 #' 'illparms.mmkin' or 'illparms.mhmkin'.
@@ -92,7 +95,7 @@ print.illparms.mmkin <- function(x, ...) {
 
 #' @rdname illparms
 #' @export
-illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, ...) {
+illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod = TRUE, slopes = TRUE, ...) {
   if (inherits(object$so, "try-error")) {
     ill_parms <- NA
   } else {
@@ -106,6 +109,15 @@ illparms.saem.mmkin <- function(object, conf.level = 0.95, random = TRUE, errmod
       ill_parms_errmod <- ints$errmod[, "lower"] < 0 & ints$errmod[, "est."] > 0
       ill_parms <- c(ill_parms, names(which(ill_parms_errmod)))
     }
+    if (slopes) {
+      if (is.null(object$so)) stop("Slope testing is only implemented for the saemix backend")
+      slope_names <- grep("^beta_", object$so@model@name.fixed, value = TRUE)
+      ci <- object$so@results@conf.int
+      rownames(ci) <- ci$name
+      slope_ci <- ci[slope_names, ]
+      ill_parms_slopes <- slope_ci[, "lower"] < 0 & slope_ci[, "estimate"] > 0
+      ill_parms <- c(ill_parms, slope_names[ill_parms_slopes])
+    }
   }
   class(ill_parms) <- "illparms.saem.mmkin"
   return(ill_parms)
diff --git a/R/intervals.R b/R/intervals.R
index 705ef6eb..fcdbaea9 100644
--- a/R/intervals.R
+++ b/R/intervals.R
@@ -78,8 +78,12 @@ intervals.saem.mmkin <- function(object, level = 0.95, backtransform = TRUE, ...
 
   # Random effects
   sdnames <- intersect(rownames(conf.int), paste("SD", pnames, sep = "."))
-  ranef_ret <- as.matrix(conf.int[sdnames, c("lower", "est.", "upper")])
-  rownames(ranef_ret) <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+  corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+  ranef_ret <- as.matrix(conf.int[c(sdnames, corrnames), c("lower", "est.", "upper")])
+  sdnames_ret <- paste0(gsub("SD\\.", "sd(", sdnames), ")")
+  corrnames_ret <- gsub("Corr\\.(.*)\\.(.*)", "corr(\\1,\\2)", corrnames)
+  rownames(ranef_ret) <- c(sdnames_ret, corrnames_ret)
+
   attr(ranef_ret, "label") <- "Random effects:"
 
 
diff --git a/R/logLik.mkinfit.R b/R/logLik.mkinfit.R
index 7cc10234..abf8517c 100644
--- a/R/logLik.mkinfit.R
+++ b/R/logLik.mkinfit.R
@@ -37,6 +37,7 @@ logLik.mkinfit <- function(object, ...) {
   val <- object$logLik
   # Number of estimated parameters
   attr(val, "df") <- length(object$bparms.optim) + length(object$errparms)
+  attr(val, "nobs") <- nobs(object)
   class(val) <- "logLik"
   return(val)
 }
diff --git a/R/mhmkin.R b/R/mhmkin.R
index 1f29dc40..6265a59e 100644
--- a/R/mhmkin.R
+++ b/R/mhmkin.R
@@ -12,13 +12,14 @@
 #' degradation models to the same data
 #' @param backend The backend to be used for fitting. Currently, only saemix is
 #' supported
-#' @param no_random_effect Default is NULL and will be passed to [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.
-#' @param auto_ranef_threshold See 'no_random_effect.
+#' @param no_random_effect Default is NULL and will be passed to [saem]. If a
+#' character vector is supplied, it will be passed to all calls to [saem],
+#' which will exclude random effects for all matching parameters. Alternatively,
+#' a list of character vectors or an object of class [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).
 #' @param algorithm The algorithm to be used for fitting (currently not used)
 #' @param \dots Further arguments that will be passed to the nonlinear mixed-effects
 #' model fitting function.
@@ -50,8 +51,44 @@ mhmkin.mmkin <- function(objects, ...) {
 
 #' @export
 #' @rdname mhmkin
+#' @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)
+#' }
 mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
-  no_random_effect = NULL, auto_ranef_threshold = 3,
+  no_random_effect = NULL,
   ...,
   cores = if (Sys.info()["sysname"] == "Windows") 1 else parallel::detectCores(), cluster = NULL)
 {
@@ -97,27 +134,42 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
   dimnames(fit_indices) <- list(degradation = names(deg_models),
                                 error = error_models)
 
+  if (is.null(no_random_effect) || is.null(dim(no_random_effect))) {
+    no_ranef <- rep(list(no_random_effect), n.fits)
+    dim(no_ranef) <- dim(fit_indices)
+  } else {
+    if (!identical(dim(no_random_effect), dim(fit_indices))) {
+      stop("Dimensions of argument 'no_random_effect' are not suitable")
+    }
+    if (is(no_random_effect, "illparms.mhmkin")) {
+      no_ranef_dim <- dim(no_random_effect)
+      no_ranef <- lapply(no_random_effect, function(x) {
+        no_ranef_split <- strsplit(x, ", ")
+        ret <- sapply(no_ranef_split, function(y) {
+          gsub("sd\\((.*)\\)", "\\1", y)
+        })
+        return(ret)
+      })
+      dim(no_ranef) <- no_ranef_dim
+    } else {
+      no_ranef <- no_random_effect
+    }
+  }
+
   fit_function <- function(fit_index) {
     w <- which(fit_indices == fit_index, arr.ind = TRUE)
     deg_index <- w[1]
     error_index <- w[2]
     mmkin_row <- objects[[error_index]][deg_index, ]
-    if (identical(no_random_effect, "auto")) {
-      ip <- illparms(mmkin_row)
-      ipt <- table(unlist(lapply(as.list(ip), strsplit, ", ")))
-      no_ranef <- names(ipt)[which((length(ds) - ipt) <= auto_ranef_threshold)]
-      transparms <- rownames(mmkin_row[[1]]$start_transformed)
-      bparms <- rownames(mmkin_row[[1]]$start)
-      names(transparms) <- bparms
-      no_ranef_trans <- transparms[no_ranef]
-    } else {
-      no_ranef_trans <- no_random_effect
-    }
     res <- try(do.call(backend_function,
-        args = c(list(mmkin_row), dot_args, list(no_random_effect = no_ranef_trans))))
+        args = c(
+          list(mmkin_row),
+          dot_args,
+          list(no_random_effect = no_ranef[[deg_index, error_index]]))))
     return(res)
   }
 
+
   fit_time <- system.time({
     if (is.null(cluster)) {
       results <- parallel::mclapply(as.list(1:n.fits), fit_function,
@@ -143,15 +195,16 @@ mhmkin.list <- function(objects, backend = "saemix", algorithm = "saem",
 #' @param j Column index selecting the fits to specific datasets
 #' @param drop If FALSE, the method always returns an mhmkin object, otherwise
 #'   either a list of fit objects or a single fit object.
-#' @return An object of class \code{\link{mhmkin}}.
+#' @return An object inheriting from \code{\link{mhmkin}}.
 #' @rdname mhmkin
 #' @export
 `[.mhmkin` <- function(x, i, j, ..., drop = FALSE) {
+  original_class <- class(x)
   class(x) <- NULL
   x_sub <- x[i, j, drop = drop]
 
   if (!drop) {
-    class(x_sub) <- "mhmkin"
+    class(x_sub) <- original_class
   }
   return(x_sub)
 }
diff --git a/R/mkinmod.R b/R/mkinmod.R
index b1fb57cb..215dbed6 100644
--- a/R/mkinmod.R
+++ b/R/mkinmod.R
@@ -26,7 +26,7 @@
 #'   Additionally, [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 [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/R/multistart.R b/R/multistart.R
index 61ef43dc..bdfbfe63 100644
--- a/R/multistart.R
+++ b/R/multistart.R
@@ -100,9 +100,12 @@ multistart.saem.mmkin <- function(object, n = 50, cores = 1,
   }
 
   if (is.null(cluster)) {
-    res <- parallel::mclapply(1:n, fit_function, mc.cores = cores)
+    res <- parallel::mclapply(1:n, fit_function,
+      mc.cores = cores, mc.preschedule = FALSE)
   } else {
-    res <- parallel::parLapply(cluster, 1:n, fit_function)
+    res <- parallel::parLapplyLB(cluster, 1:n, fit_function,
+      mc.preschedule = FALSE
+    )
   }
   attr(res, "orig") <- object
   attr(res, "start_parms") <- start_parms
diff --git a/R/parms.R b/R/parms.R
index bd4e479b..bb04a570 100644
--- a/R/parms.R
+++ b/R/parms.R
@@ -77,6 +77,6 @@ parms.multistart <- function(object, exclude_failed = TRUE, ...) {
   successful <- which(!is.na(res[, 1]))
   first_success <- successful[1]
   colnames(res) <- names(parms(object[[first_success]]))
-  if (exclude_failed) res <- res[successful, ]
+  if (exclude_failed[1]) res <- res[successful, ]
   return(res)
 }
diff --git a/R/parplot.R b/R/parplot.R
index 627a4eb8..3da4b51a 100644
--- a/R/parplot.R
+++ b/R/parplot.R
@@ -4,9 +4,16 @@
 #' 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).
 #'
+#' 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.
+#'
 #' @param object The [multistart] object
 #' @param llmin The minimum likelihood of objects to be shown
-#' @param scale By default, scale parameters using the best available fit.
+#' @param llquant Fractional value for selecting only the fits with higher
+#' likelihoods. Overrides 'llmin'.
+#' @param scale By default, scale parameters using the best
+#' available fit.
 #' If 'median', parameters are scaled using the median parameters from all fits.
 #' @param main Title of the plot
 #' @param lpos Positioning of the legend.
@@ -16,7 +23,7 @@
 #' of the in vitro erythropoiesis. BMC Bioinformatics. 2021 Oct 4;22(1):478.
 #' doi: 10.1186/s12859-021-04373-4.
 #' @seealso [multistart]
-#' @importFrom stats median
+#' @importFrom stats median quantile
 #' @export
 parplot <- function(object, ...) {
   UseMethod("parplot")
@@ -24,7 +31,8 @@ parplot <- function(object, ...) {
 
 #' @rdname parplot
 #' @export
-parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best", "median"),
+parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, llquant = NA,
+  scale = c("best", "median"),
   lpos = "bottomleft", main = "", ...)
 {
   oldpar <- par(no.readonly = TRUE)
@@ -32,8 +40,8 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
 
   orig <- attr(object, "orig")
   orig_parms <- parms(orig)
-  start_parms <- orig$mean_dp_start
-  all_parms <- parms(object)
+  start_degparms <- orig$mean_dp_start
+  all_parms <- parms(object, exclude_failed = FALSE)
 
   if (inherits(object, "multistart.saem.mmkin")) {
     llfunc <- function(object) {
@@ -44,23 +52,27 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
     stop("parplot is only implemented for multistart.saem.mmkin objects")
   }
   ll <- sapply(object, llfunc)
+  if (!is.na(llquant[1])) {
+    if (llmin != -Inf) warning("Overriding 'llmin' because 'llquant' was specified")
+    llmin <- quantile(ll, 1 - llquant)
+  }
   selected <- which(ll > llmin)
   selected_parms <- all_parms[selected, ]
 
   par(las = 1)
   if (orig$transformations == "mkin") {
-    degparm_names_transformed <- names(start_parms)
+    degparm_names_transformed <- names(start_degparms)
     degparm_index <- which(names(orig_parms) %in% degparm_names_transformed)
     orig_parms[degparm_names_transformed] <- backtransform_odeparms(
       orig_parms[degparm_names_transformed],
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
-    start_parms <- backtransform_odeparms(start_parms,
+    start_degparms <- backtransform_odeparms(start_degparms,
       orig$mmkin[[1]]$mkinmod,
       transform_rates = orig$mmkin[[1]]$transform_rates,
       transform_fractions = orig$mmkin[[1]]$transform_fractions)
-    degparm_names <- names(start_parms)
+    degparm_names <- names(start_degparms)
 
     names(orig_parms) <- c(degparm_names, names(orig_parms[-degparm_index]))
 
@@ -72,6 +84,13 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
     colnames(selected_parms)[1:length(degparm_names)] <- degparm_names
   }
 
+  start_errparms <- orig$so@model@error.init
+  names(start_errparms) <- orig$so@model@name.sigma
+
+  start_omegaparms <- orig$so@model@omega.init
+
+  start_parms <- c(start_degparms, start_errparms)
+
   scale <- match.arg(scale)
   parm_scale <- switch(scale,
     best = selected_parms[which.best(object[selected]), ],
@@ -99,7 +118,7 @@ parplot.multistart.saem.mmkin <- function(object, llmin = -Inf, scale = c("best"
   legend(lpos, inset = c(0.05, 0.05), bty = "n",
     pch = 1, col = 3:1, lty = c(NA, NA, 1),
     legend = c(
-      "Starting parameters",
-      "Original run",
+      "Original start",
+      "Original results",
       "Multistart runs"))
 }
diff --git a/R/read_spreadsheet.R b/R/read_spreadsheet.R
index a20af6db..7ad09c3e 100644
--- a/R/read_spreadsheet.R
+++ b/R/read_spreadsheet.R
@@ -37,7 +37,7 @@
 #' and moisture normalisation factors in the sheet 'Datasets'?
 #' @export
 read_spreadsheet <- function(path, valid_datasets = "all",
-  parent_only = TRUE, normalize = TRUE)
+  parent_only = FALSE, normalize = TRUE)
 {
   if (!requireNamespace("readxl", quietly = TRUE))
     stop("Please install the readxl package to use this function")
diff --git a/R/saem.R b/R/saem.R
index 5b8021de..b29cf8a9 100644
--- a/R/saem.R
+++ b/R/saem.R
@@ -149,7 +149,7 @@ saem.mmkin <- function(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,
@@ -718,6 +718,7 @@ saemix_model <- function(object, solution_type = "auto",
     covariance.model = covariance.model,
     covariate.model = covariate.model,
     omega.init = omega.init,
+    error.init = error.init,
     ...
   )
   attr(res, "mean_dp_start") <- degparms_optim
diff --git a/R/summary.saem.mmkin.R b/R/summary.saem.mmkin.R
index 2754e9f0..49b02a50 100644
--- a/R/summary.saem.mmkin.R
+++ b/R/summary.saem.mmkin.R
@@ -75,10 +75,21 @@
 #' 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)
 #' }
 #'
 #' @export
@@ -136,10 +147,11 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
   }
 
   #  Correlation of fixed effects (inspired by summary.nlme)
-  varFix <- try(vcov(object$so)[1:n_fixed, 1:n_fixed])
-  if (inherits(varFix, "try-error")) {
+  cov_so <- try(solve(object$so@results@fim), silent = TRUE)
+  if (inherits(cov_so, "try-error")) {
     object$corFixed <- NA
   } else {
+    varFix <- cov_so[1:n_fixed, 1:n_fixed]
     stdFix <- sqrt(diag(varFix))
     object$corFixed <- array(
       t(varFix/stdFix)/stdFix,
@@ -149,7 +161,8 @@ summary.saem.mmkin <- function(object, data = FALSE, verbose = FALSE, distimes =
 
   # Random effects
   sdnames <- intersect(rownames(conf.int), paste0("SD.", pnames))
-  confint_ranef <- as.matrix(conf.int[sdnames, c("estimate", "lower", "upper")])
+  corrnames <- grep("^Corr.", rownames(conf.int), value = TRUE)
+  confint_ranef <- as.matrix(conf.int[c(sdnames, corrnames), c("estimate", "lower", "upper")])
   colnames(confint_ranef)[1] <- "est."
 
   # Error model
@@ -225,13 +238,22 @@ print.summary.saem.mmkin <- function(x, digits = max(3, getOption("digits") - 3)
     obs = "Variance unique to each observed variable",
     tc = "Two-component variance function"), "\n")
 
-  cat("\nMean of starting values for individual parameters:\n")
+  cat("\nStarting values for degradation parameters:\n")
   print(x$mean_dp_start, digits = digits)
 
   cat("\nFixed degradation parameter values:\n")
   if(length(x$fixed$value) == 0) cat("None\n")
   else print(x$fixed, digits = digits)
 
+  cat("\nStarting values for random effects (square root of initial entries in omega):\n")
+  print(sqrt(x$so@model@omega.init), digits = digits)
+
+  cat("\nStarting values for error model parameters:\n")
+  errparms <- x$so@model@error.init
+  names(errparms) <- x$so@model@name.sigma
+  errparms <- errparms[x$so@model@indx.res]
+  print(errparms, digits = digits)
+
   cat("\nResults:\n\n")
   cat("Likelihood computed by importance sampling\n")
   print(data.frame(AIC = x$AIC, BIC = x$BIC, logLik = x$logLik,
diff --git a/R/summary_listing.R b/R/summary_listing.R
new file mode 100644
index 00000000..38b52240
--- /dev/null
+++ b/R/summary_listing.R
@@ -0,0 +1,59 @@
+#' Display the output of a summary function according to the output format
+#'
+#' This function is intended for use in a R markdown code chunk with the chunk
+#' option `results = "asis"`.
+#'
+#' @param object The object for which the summary is to be listed
+#' @param caption An optional caption
+#' @param label An optional label, ignored in html output
+#' @param clearpage Should a new page be started after the listing? Ignored in html output
+#' @export
+summary_listing <- function(object, caption = NULL, label = NULL,
+  clearpage = TRUE) {
+  if (knitr::is_latex_output()) {
+    tex_listing(object = object, caption = caption, label = label,
+      clearpage = clearpage)
+  }
+  if (knitr::is_html_output()) {
+    html_listing(object = object, caption = caption)
+  }
+}
+
+#' @rdname summary_listing
+#' @export
+tex_listing <- function(object, caption = NULL, label = NULL,
+  clearpage = TRUE) {
+  cat("\n")
+  cat("\\begin{listing}", "\n")
+  if (!is.null(caption)) {
+    cat("\\caption{", caption, "}", "\n", sep = "")
+  }
+  if (!is.null(label)) {
+    cat("\\caption{", label, "}", "\n", sep = "")
+  }
+  cat("\\begin{snugshade}", "\n")
+  cat("\\scriptsize", "\n")
+  cat("\\begin{verbatim}", "\n")
+  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
+  cat("\n")
+  cat("\\end{verbatim}", "\n")
+  cat("\\end{snugshade}", "\n")
+  cat("\\end{listing}", "\n")
+  if (clearpage) {
+    cat("\\clearpage", "\n")
+  }
+}
+
+#' @rdname summary_listing
+#' @export
+html_listing <- function(object, caption = NULL) {
+  cat("\n")
+  if (!is.null(caption)) {
+    cat("<caption>", caption, "</caption>", "\n", sep = "")
+  }
+  cat("<pre><code>\n")
+  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
+  cat("\n")
+  cat("</pre></code>\n")
+}
+
diff --git a/R/tex_listing.R b/R/tex_listing.R
deleted file mode 100644
index 05f662e4..00000000
--- a/R/tex_listing.R
+++ /dev/null
@@ -1,32 +0,0 @@
-#' Wrap the output of a summary function in tex listing environment
-#'
-#' This function can be used in a R markdown code chunk with the chunk
-#' option `results = "asis"`.
-#'
-#' @param object The object for which the summary is to be listed
-#' @param caption An optional caption
-#' @param label An optional label
-#' @param clearpage Should a new page be started after the listing?
-#' @export
-tex_listing <- function(object, caption = NULL, label = NULL,
-  clearpage = TRUE) {
-  cat("\n")
-  cat("\\begin{listing}", "\n")
-  if (!is.null(caption)) {
-    cat("\\caption{", caption, "}", "\n", sep = "")
-  }
-  if (!is.null(label)) {
-    cat("\\caption{", label, "}", "\n", sep = "")
-  }
-  cat("\\begin{snugshade}", "\n")
-  cat("\\scriptsize", "\n")
-  cat("\\begin{verbatim}", "\n")
-  cat(capture.output(suppressWarnings(summary(object))), sep = "\n")
-  cat("\n")
-  cat("\\end{verbatim}", "\n")
-  cat("\\end{snugshade}", "\n")
-  cat("\\end{listing}", "\n")
-  if (clearpage) {
-    cat("\\clearpage", "\n")
-  }
-}