diff options
Diffstat (limited to 'R')
| -rw-r--r-- | R/ds_mixed.R | 17 | ||||
| -rw-r--r-- | R/hierarchical_kinetics.R | 40 | ||||
| -rw-r--r-- | R/illparms.R | 14 | ||||
| -rw-r--r-- | R/intervals.R | 8 | ||||
| -rw-r--r-- | R/logLik.mkinfit.R | 1 | ||||
| -rw-r--r-- | R/mhmkin.R | 97 | ||||
| -rw-r--r-- | R/mkinmod.R | 2 | ||||
| -rw-r--r-- | R/multistart.R | 7 | ||||
| -rw-r--r-- | R/parms.R | 2 | ||||
| -rw-r--r-- | R/parplot.R | 39 | ||||
| -rw-r--r-- | R/read_spreadsheet.R | 2 | ||||
| -rw-r--r-- | R/saem.R | 3 | ||||
| -rw-r--r-- | R/summary.saem.mmkin.R | 32 | ||||
| -rw-r--r-- | R/summary_listing.R | 59 | ||||
| -rw-r--r-- | R/tex_listing.R | 32 |
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) } @@ -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 @@ -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") @@ -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") - } -} |