diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml index 27e1293..d0f6bbd 100644 --- a/.github/workflows/R-CMD-check.yaml +++ b/.github/workflows/R-CMD-check.yaml @@ -6,7 +6,9 @@ on: pull_request: branches: [main, master] -name: R-CMD-check +name: R-CMD-check.yaml + +permissions: read-all jobs: R-CMD-check: @@ -22,17 +24,14 @@ jobs: - {os: windows-latest, r: 'release'} - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} - {os: ubuntu-latest, r: 'release'} - #- {os: ubuntu-latest, r: 'oldrel-1'} + - {os: ubuntu-latest, r: 'oldrel-1'} env: - R_REMOTES_NO_ERRORS_FROM_WARNINGS: true - RSPM: ${{ matrix.config.rspm }} GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - _R_CHECK_FORCE_SUGGESTS_: false R_KEEP_PKG_SOURCE: yes steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 @@ -44,7 +43,6 @@ jobs: - uses: r-lib/actions/setup-r-dependencies@v2 with: - pak-version: devel extra-packages: | any::rcmdcheck nlmixr2/lotri @@ -52,8 +50,13 @@ jobs: nlmixr2/rxode2 nlmixr2/nlmixr2est nlmixr2/nlmixr2plot + nlmixr2/PreciseSums + nlmixr2/lbfgsb3c + nlmixr2/n1qn1c + nlmixr2/dparser-R needs: check - uses: r-lib/actions/check-r-package@v2 with: upload-snapshots: true + build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index f8afb66..b4bbaf3 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -9,7 +9,9 @@ on: types: [published] workflow_dispatch: -name: pkgdown +name: pkgdown.yaml + +permissions: read-all jobs: pkgdown: @@ -19,8 +21,10 @@ jobs: group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} env: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} + permissions: + contents: write steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-pandoc@v2 @@ -37,6 +41,10 @@ jobs: nlmixr2/rxode2 nlmixr2/nlmixr2est nlmixr2/nlmixr2plot + nlmixr2/PreciseSums + nlmixr2/lbfgsb3c + nlmixr2/n1qn1c + nlmixr2/dparser-R local::. needs: website @@ -46,7 +54,7 @@ jobs: - name: Deploy to GitHub pages 🚀 if: github.event_name != 'pull_request' - uses: JamesIves/github-pages-deploy-action@v4.4.1 + uses: JamesIves/github-pages-deploy-action@v4.5.0 with: clean: false branch: gh-pages diff --git a/.github/workflows/test-coverage.yaml b/.github/workflows/test-coverage.yaml index 4907ffa..ca6b300 100644 --- a/.github/workflows/test-coverage.yaml +++ b/.github/workflows/test-coverage.yaml @@ -6,7 +6,9 @@ on: pull_request: branches: [main, master] -name: test-coverage +name: test-coverage.yaml + +permissions: read-all jobs: test-coverage: @@ -15,7 +17,7 @@ jobs: GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} steps: - - uses: actions/checkout@v3 + - uses: actions/checkout@v4 - uses: r-lib/actions/setup-r@v2 with: @@ -25,32 +27,46 @@ jobs: with: extra-packages: | any::covr + any::xml2 nlmixr2/lotri nlmixr2/rxode2ll nlmixr2/rxode2 nlmixr2/nlmixr2est nlmixr2/nlmixr2plot + nlmixr2/PreciseSums + nlmixr2/lbfgsb3c + nlmixr2/n1qn1c + nlmixr2/dparser-R needs: coverage - name: Test coverage run: | - covr::codecov( + cov <- covr::package_coverage( quiet = FALSE, clean = FALSE, - install_path = file.path(Sys.getenv("RUNNER_TEMP"), "package") + install_path = file.path(normalizePath(Sys.getenv("RUNNER_TEMP"), winslash = "/"), "package") ) + covr::to_cobertura(cov) shell: Rscript {0} + - uses: codecov/codecov-action@v4 + with: + fail_ci_if_error: ${{ github.event_name != 'pull_request' && true || false }} + file: ./cobertura.xml + plugin: noop + disable_search: true + token: ${{ secrets.CODECOV_TOKEN }} + - name: Show testthat output if: always() run: | ## -------------------------------------------------------------------- - find ${{ runner.temp }}/package -name 'testthat.Rout*' -exec cat '{}' \; || true + find '${{ runner.temp }}/package' -name 'testthat.Rout*' -exec cat '{}' \; || true shell: bash - name: Upload test results if: failure() - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: coverage-test-failures path: ${{ runner.temp }}/package diff --git a/DESCRIPTION b/DESCRIPTION index 74fe1e0..535c3c0 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: nlmixr2extra Title: Nonlinear Mixed Effects Models in Population PK/PD, Extra Support Functions -Version: 2.0.10 +Version: 3.0.0 Authors@R: c( person("Matthew", "Fidler", email="matthew.fidler@gmail.com", role = c("aut", "cre"), comment = c(ORCID = "0000-0001-8538-6691")), @@ -32,6 +32,7 @@ Imports: crayon, data.table, digest, + dplyr, ggplot2, ggtext, knitr, @@ -49,7 +50,6 @@ Suggests: nlmixr2data, testthat (>= 3.0.0), withr, - dplyr, devtools LinkingTo: Rcpp, @@ -59,5 +59,5 @@ Config/testthat/edition: 3 Encoding: UTF-8 NeedsCompilation: yes Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 LazyData: true diff --git a/NAMESPACE b/NAMESPACE index c7b34a7..f6e57b5 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,9 +1,21 @@ # Generated by roxygen2: do not edit by hand S3method(bootplot,nlmixr2FitCore) +S3method(extractEqHelper,"(") +S3method(extractEqHelper,"function") +S3method(extractEqHelper,"{") +S3method(extractEqHelper,call) +S3method(extractEqHelper,character) +S3method(extractEqHelper,default) +S3method(extractEqHelper,name) +S3method(extractEqHelper,nlmixr2FitCore) +S3method(extractEqHelper,numeric) +S3method(extractEqHelper,rxUi) S3method(knit_print,nlmixr2FitCore) S3method(knit_print,rxUi) S3method(print,nlmixr2BoostrapSummary) +S3method(profile,nlmixr2FitCore) +S3method(rxUiDeparse,llpControl) export(adaptivelassoCoefficients) export(addCatCovariates) export(addorremoveCovariate) @@ -14,6 +26,7 @@ export(bootstrapFit) export(buildcovInfo) export(buildupatedUI) export(covarSearchAuto) +export(fixedControl) export(foldgen) export(forwardSearch) export(horseshoeSummardf) @@ -21,6 +34,7 @@ export(ini) export(knit_print) export(lassoCoefficients) export(lassoSummardf) +export(llpControl) export(model) export(nlmixr) export(nlmixr2) @@ -28,7 +42,11 @@ export(nlmixrWithTiming) export(normalizedData) export(optimUnisampling) export(preconditionFit) +export(profileFixed) +export(profileFixedSingle) +export(profileLlp) export(regularmodel) +export(rxUiDeparse) import(lotri) import(utils) importFrom(Rcpp,evalCpp) @@ -40,6 +58,7 @@ importFrom(nlmixr2est,nlmixrWithTiming) importFrom(nlmixr2est,setCov) importFrom(rxode2,ini) importFrom(rxode2,model) +importFrom(rxode2,rxUiDeparse) importFrom(stats,AIC) importFrom(stats,approxfun) importFrom(stats,cov) diff --git a/NEWS.md b/NEWS.md index 6f26dab..eaace35 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# nlmixr2extra development version + +* New `profile()` method for likelihood profiling (Issue #1) + # nlmixr2extra 2.0.10 * `bootstrapFit()` fixes `se` option (Issue #66) diff --git a/R/computingutil.R b/R/computingutil.R index 783318c..995eae5 100644 --- a/R/computingutil.R +++ b/R/computingutil.R @@ -958,11 +958,7 @@ modelBootstrap <- function(fit, } # get control settings for the 'fit' object and save computation effort by not computing the tables - .ctl <- fit$control - .ctl$print <- 0L - .ctl$covMethod <- 0L - .ctl$calcTables <- FALSE - .ctl$compress <- FALSE + .ctl <- setQuietFastControl(fit$control) modelsEnsemble <- lapply(bootData[.env$mod_idx:nboot], function(boot_data) { diff --git a/R/knit_printEquation.R b/R/knit_printEquation.R index 696d41e..faf8c28 100644 --- a/R/knit_printEquation.R +++ b/R/knit_printEquation.R @@ -69,22 +69,27 @@ extractEqHelperLhsRhs <- function(x, ..., inModel, prefix = "", middle, suffix = paste0(prefix, lhs, middle, rhs, suffix) } +#' @export extractEqHelper.nlmixr2FitCore <- function(x, ..., inModel) { extractEqHelper(as.function(x), ..., inModel=FALSE) } +#' @export extractEqHelper.rxUi <- function(x, ..., inModel, name) { extractEqHelper(as.function(x), ..., inModel=FALSE) } +#' @export extractEqHelper.function <- function(x, ..., inModel, name) { extractEqHelper(methods::functionBody(x), ..., inModel=FALSE) } +#' @export "extractEqHelper.{" <- function(x, ..., inModel, name) { extractEqHelperSeqDrop(x, ..., inModel = inModel) } +#' @export "extractEqHelper.(" <- function(x, ..., inModel, dropParen = FALSE, name) { stopifnot(length(x) == 2) if (dropParen) { @@ -178,6 +183,7 @@ latexOpMap <- "!"="{\\lnot}" ) +#' @export extractEqHelper.call <- function(x, ..., inModel, name) { if (inModel) { if (is.name(x[[1]])) { @@ -270,6 +276,7 @@ extractEqHelper.call <- function(x, ..., inModel, name) { ret } +#' @export extractEqHelper.name <- function(x, ..., inModel, underscoreToSubscript = FALSE, name) { if (inModel) { ret <- as.character(x) @@ -292,6 +299,7 @@ extractEqHelper.name <- function(x, ..., inModel, underscoreToSubscript = FALSE, ret } +#' @export extractEqHelper.numeric <- function(x, ..., inModel, name = NULL) { if (inModel) { ret <- format(x) @@ -329,6 +337,7 @@ escapeLatex <- function (x, newlines = FALSE, spaces = FALSE) { x } +#' @export extractEqHelper.character <- function(x, ..., inModel, name = NULL) { if (inModel) { ret <- sprintf('\\text{"%s"}', escapeLatex(x)) @@ -394,6 +403,7 @@ extractEqHelper.if <- function(x, ..., inModel, alignment, indent = 0L, firstIf ret } +#' @export extractEqHelper.default <- function(x, ..., inModel) { if (inherits(x, "<-") | inherits(x, "=")) { # The assignment classes go via extractEqHelper.default to fix an R CMD diff --git a/R/profile.R b/R/profile.R new file mode 100644 index 0000000..447db1f --- /dev/null +++ b/R/profile.R @@ -0,0 +1,463 @@ +# General profiling functions ---- + +#' Perform likelihood profiling on nlmixr2 focei fits +#' +#' +#' @details +#' +#' # Log-likelihood profiling +#' +#' `method = "llp"` +#' +#' The search will stop when either the OFV is within `ofvtol` of the desired +#' OFV change or when the parameter is interpolating to more significant digits +#' than specified in `paramDigits`. The "llp" method uses the `profileLlp()` +#' function. See its help for more details. +#' +#' # Fixed points +#' +#' `method = "fixed"` +#' +#' Estimate the OFV for specific fixed values. The "fixed" method uses the +#' `profileFixed()` function. See its help for more details. +#' +#' @param fitted The fit model +#' @param ... ignored +#' @param which The parameter names to perform likelihood profiling on +#' (`NULL` indicates all parameters) +#' @param method Method to use for profiling (see the details) +#' @param control Control arguments for the `method` +#' @return A data.frame with one column named `Parameter` indicating the +#' parameter being fixed on that row, one column for the `OFV` indicating the +#' OFV estimated for the model at that step, one column named `profileBound` +#' indicating the estimated value for the profile likelihood and its step +#' above the minimum profile likelihood value, and columns for each parameter +#' estimate (or fixed) in the model. +#' @family Profiling +#' @examples +#' \dontrun{ +#' # Likelihood profiling takes a long time to run each model multiple times, so +#' # be aware that running this example may take a few minutes. +#' oneCmt <- function() { +#' ini({ +#' tka <- log(1.57) +#' tcl <- log(2.72) +#' tv <- fixed(log(31.5)) +#' eta.ka ~ 0.6 +#' add.sd <- 0.7 +#' }) +#' model({ +#' ka <- exp(tka + eta.ka) +#' cl <- exp(tcl) +#' v <- exp(tv) +#' cp <- linCmt() +#' cp ~ add(add.sd) +#' }) +#' } +#' +#' fit <- +#' nlmixr2( +#' oneCmt, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) +#' ) +#' # profile all parameters +#' profall <- profile(fit) +#' +#' # profile a single parameter +#' proftka <- profile(fit, which = "tka") +#' } +#' @export +profile.nlmixr2FitCore <- function(fitted, ..., + which = NULL, + method = c("llp", "fixed"), + control = list()) { + method <- match.arg(method) + + if (method == "llp") { + profileLlp(fitted = fitted, which = which, control = control) + } else if (method == "fixed") { + profileFixed(fitted = fitted, which = which, control = control) + } else { + stop("Invalid 'method': ", method) # nocov + } +} + +#' Give the output data.frame for a single model for profile.nlmixr2FitCore +#' +#' @inheritParams profile.nlmixr2FitCore +#' @param fixedVal The value that `which` is fixed to in case the model does not +#' converge. +#' @return A data.frame with columns named "Parameter" (the parameter name(s) +#' that were fixed), OFV (the objective function value), and the current +#' estimate for each of the parameters. Omega values are given as their +#' variances and covariances. +#' @family Profiling +profileNlmixr2FitCoreRet <- function(fitted, which, fixedVal) { + if (inherits(fitted, "try-error")) { + ret <- data.frame(Parameter = which, OFV = NA_real_, X = fixedVal) + names(ret)[3] <- which + } else { + if (any(names(nlmixr2est::fixef(fitted)) %in% c("Parameter", "OFV"))) { + cli::cli_abort("Cannot profile a model with a parameter named either 'Parameter' or 'OFV'") + } + retOmega <- as.data.frame(t(diag(fitted$omega))) + for (idxRow in seq_len(nrow(fitted$omega)) - 1) { + for (idxCol in seq_len(idxRow)) { + # Only capture the columns when the covariance is nonzero + if (fitted$omega[idxRow + 1, idxCol] != 0) { + covColName <- + sprintf( + "cov(%s,%s)", + rownames(fitted$omega)[idxRow + 1], + colnames(fitted$omega)[idxCol] + ) + retOmega[[covColName]] <- fitted$omega[idxRow + 1, idxCol] + } + } + } + ret <- + cbind( + data.frame(Parameter = which, OFV = fitted$objective), + data.frame(t(nlmixr2est::fixef(fitted))), + retOmega + ) + } + rownames(ret) <- NULL + ret +} + +# Fixed parameter estimate profiling ---- + +#' Estimate the objective function values for a model while fixing defined +#' parameter values +#' +#' @inheritParams profile.nlmixr2FitCore +#' @param which A data.frame with column names of parameters to fix and values +#' of the fitted value to fix (one row per set of parameters to estimate) +#' @param control A list passed to `fixedControl()` (currently unused) +#' @inherit profileNlmixr2FitCoreRet return +#' @family Profiling +#' @export +profileFixed <- function(fitted, which, control = list()) { + control <- do.call(fixedControl, control) + checkmate::assert_data_frame(which, types = "numeric", any.missing = FALSE, min.rows = 1) + dplyr::bind_rows(lapply( + X = seq_len(nrow(which)), + FUN = \(idx, fitted) profileFixedSingle(fitted = fitted, which = which[idx, , drop = FALSE]), + fitted = fitted + )) +} + +#' @describeIn profileFixed Estimate the objective function value for a model +#' while fixing a single set of defined parameter values (for use in parameter +#' profiling) +#' +#' @param which A data.frame with column names of parameters to fix and values +#' of the fitted value to fix (one row only). +#' @returns `which` with a column named `OFV` added with the objective function +#' value of the fitted estimate fixing the parameters in the other columns +#' @family Profiling +#' @export +profileFixedSingle <- function(fitted, which) { + checkmate::assert_data_frame(which, types = "numeric", any.missing = FALSE, nrows = 1) + checkmate::assert_names(names(which), subset.of = names(nlmixr2est::fixef(fitted))) + + message("Profiling ", paste(names(which), "=", unlist(which), collapse = ", ")) + # Update the model by fixing all of the parameters + paramToFix <- lapply(X = which, FUN = \(x) str2lang(sprintf("fixed(%s)", x))) + iniArgs <- append(list(x = fitted), paramToFix) + modelToFit <- suppressMessages(do.call(rxode2::ini, iniArgs)) + controlFit <- fitted$control + newFit <- + try(suppressMessages( + nlmixr2est::nlmixr2( + modelToFit, + est = fitted$est, + control = setQuietFastControl(fitted$control) + ) + )) + ret <- profileNlmixr2FitCoreRet(fitted = newFit, which = paste(names(which), collapse = ",")) + message("OFV: ", ret$OFV) + ret +} + +#' Control options for fixed-value likelihood profiling +#' +#' @returns A validated list of control options for fixed-value likelihood +#' profiling +#' @family Profiling +#' @seealso [profileFixed()] +#' @export +fixedControl <- function() { + ret <- list() + class(ret) <- "fixedControl" + ret +} + +# Log-likelihood profiling ---- + +#' Profile confidence intervals with log-likelihood profiling +#' +#' @inheritParams profile.nlmixr2FitCore +#' @param control A list passed to `llpControl()` +#' @param which Either `NULL` to profile all parameters or a character vector of +#' parameters to estimate +#' @return A data.frame with columns named "Parameter" (the parameter name(s) +#' that were fixed), OFV (the objective function value), and the current +#' estimate for each of the parameters. In addition, if any boundary is +#' found, the OFV increase will be indicated by the absolute value of the +#' "profileBound" column and if that boundary is the upper or lower boundary +#' will be indicated by the "profileBound" column being positive or negative, +#' respectively. +#' @family Profiling +#' @export +profileLlp <- function(fitted, which, control) { + # Validate inputs + control <- do.call(llpControl, control) + + if (is.null(which)) { + which <- names(nlmixr2est::fixef(fitted)) + } else { + checkmate::assert_subset(x = which, choices = names(nlmixr2est::fixef(fitted)), empty.ok = FALSE) + } + + if (length(which) > 1) { + ret <- data.frame() + # Make the general case of profiling many parameters to be a concatenation + # of single cases + for (currentWhich in which) { + ret <- + dplyr::bind_rows( + ret, + profileLlp( + fitted = fitted, + which = currentWhich, + control = control + ) + ) + } + } else { + if (length(control$rseTheta) == 1 & is.null(names(control$rseTheta))) { + control$rseTheta <- setNames(rep(control$rseTheta, length(which)), which) + } else { + checkmate::assert_names(names(control$rseTheta), subset.of = names(nlmixr2est::fixef(fitted))) + checkmate::assert_numeric(control$rseTheta, lower = 0, any.missing = FALSE, min.len = 1) + } + + # The original fitted model should be in the output + ret <- profileNlmixr2FitCoreRet(fitted = fitted, which = which) + + effectRange <- fitted$iniDf[fitted$iniDf$name == which, c("lower", "upper")] + estInitial <- + profileNlmixr2FitDataEstInitial( + estimates = ret, + which = which, + ofvIncrease = control$ofvIncrease, + rseTheta = control$rseTheta, + lower = effectRange$lower, + upper = effectRange$upper + ) + # Find the lower and upper limits + for (direction in c(-1, 1)) { + if (direction == -1) { + currentInitialEst = estInitial[1] + currentUpper <- nlmixr2est::fixef(fitted)[[which]] + currentLower <- effectRange$lower + } else { + currentInitialEst = estInitial[2] + currentLower <- nlmixr2est::fixef(fitted)[[which]] + currentUpper <- effectRange$upper + } + + retCurrentDirection <- + optimProfile( + par = currentInitialEst, + fitted = fitted, + which = which, + optimDf = ret, + direction = direction, + ofvIncrease = control$ofvIncrease, + lower = currentLower, upper = currentUpper, + itermax = control$itermax, + ofvtol = control$ofvtol, + paramDigits = control$paramDigits + ) + ret <- unique(dplyr::bind_rows(ret, retCurrentDirection)) + } + } + ret +} + +# Single-parameter, single-direction optimization of likelihood profiling +optimProfile <- function(par, fitted, optimDf, which, ofvIncrease, direction, lower = -Inf, upper = Inf, + itermax = 10, ofvtol = 0.005, paramDigits = 3, + extrapolateExpand = 1.5) { + currentOFVDiff <- Inf + currentIter <- 0 + currentPar <- par + origMinOFV <- min(optimDf$OFV, na.rm = TRUE) + ret <- optimDf[optimDf$Parameter %in% which & optimDf[[which]] >= lower & optimDf[[which]] <= upper, ] + while (itermax > currentIter) { + currentIter <- currentIter + 1 + dfWhich <- stats::setNames(data.frame(X = currentPar), nm = which) + fitResult <- profileFixedSingle(fitted = fitted, which = dfWhich) + + ret <- dplyr::bind_rows(ret, fitResult) + ret <- ret[order(ret$OFV), , drop = FALSE] + currentMinOFV <- min(ret$OFV, na.rm = TRUE) + if (isTRUE(any(currentMinOFV < origMinOFV))) { + warning("OFV decreased while profiling ", which, " profile may be inaccurate") + origMinOFV <- NA # Only give the warning once + } + targetOfv <- currentMinOFV + ofvIncrease + if (all(targetOfv > ret$OFV)) { # extrapolate + # Find the closest two rows for extrapolation + extrapRows <- ret[!is.na(ret$OFV), ] + if (nrow(extrapRows) < 2) { + messageProfileComplete(which, direction = direction, "aborted due to lack of convergence prior to attempted extrapolation") + currentIter <- Inf + } + # When generating the extrapolation rows, always have the first row as the + # most extreme (farthest out) + if (direction > 0) { + extrapRows <- ret[nrow(ret) + c(-1, 0), ] + } else { + extrapRows <- ret[1:2, ] + } + currentPar <- + extrapolateExpand*diff(extrapRows[[which]])/diff(extrapRows$OFV)*(targetOfv - extrapRows$OFV[1]) + extrapRows[[which]][1] + # ensure that we are within the boundaries with a slight margin + margin <- sqrt(.Machine$double.eps) + currentPar <- min(max(currentPar, lower + margin), upper - margin) + } else { # interpolate + currentPar <- stats::approx(x = ret$OFV, y = ret[[which]], xout = targetOfv)$y + } + if (currentPar %in% ret[[which]]) { + # Don't try to test the same parameter value multiple times (usually the + # case for hitting a limit) + messageProfileComplete(which, direction = direction, "aborted due to attempted repeated parameter value estimation") + currentIter <- Inf + } + # Check if the parameter significant digits are sufficiently precise. This + # is done by finding values that have the same significant digit rounding + # where one estimate is above and one is below the next parameter to be + # estimated. + signifSame <- signif(ret[[which]], digits = paramDigits) == signif(currentPar, digits = paramDigits) + signifAbove <- ret[[which]] > currentPar + signifBelow <- ret[[which]] < currentPar + convergedSignif <- any(signifSame & signifAbove) && any(signifSame & signifBelow) + + currentOFVDiff <- min(abs(targetOfv - ret$OFV), na.rm = TRUE) + convergedOfv <- currentOFVDiff <= ofvtol + + converged <- convergedSignif | convergedOfv + + if (converged) { + if (convergedSignif & convergedOfv) { + messageProfileComplete(which, direction = direction, "complete due to significant digits precision and achieving OFV within specified tolerance") + } else if (convergedSignif) { + messageProfileComplete(which, direction = direction, "complete due to significant digits precision") + } else if (convergedOfv) { + messageProfileComplete(which, direction = direction, "complete due to achieving OFV within specified tolerance") + } else { + stop("Unclear convergece criteria, please report a bug") # nocov + } + storeBound <- data.frame(Parameter = which, X = currentPar, profileBound = direction*ofvIncrease) + names(storeBound)[2] <- which + ret <- dplyr::bind_rows(storeBound, ret) + currentIter <- Inf + } else { + message("Profiling ", which, " best distance to target OFV: ", currentOFVDiff) + } + } + if (is.finite(currentIter)) { + messageProfileComplete(which, direction = direction, "aborted due to too many iterations without convergence") + } + + # Sort the output by parameter value + ret <- ret[order(ret[[which]]), ] + # Add extra information to assist with subsequent + # The output rownames are uninformative and distracting; remove them. + rownames(ret) <- NULL + ret +} + +#' Control options for log-likelihood profiling +#' +#' @param ofvIncrease The targetted change in objective function value (3.84 +#' corresponds to a Chi-squared test with a 95% confidence interval) +#' @param rseTheta The relative standard error (percent) for the model +#' parameters. It can be missing (the default) in which case a default value +#' of 30% will be applied. If given as a single number, it will be applied to +#' all parameters. If given as a named vector of numbers, it will be applied +#' to each named parameter. +#' @param itermax Maximum number of likelihood profiling iterations for each +#' bound estimated +#' @param ofvtol The relative tolerance for the objective function being close +#' enough to the `ofvIncrease`. +#' @param paramDigits The number of significant digits required for the +#' parameter. When interpolation attempts to get smaller than that number of +#' significant digits, it will stop. +#' @param extrapolateExpand When extrapolating outside the range previously +#' tested, how far should the step occur as a ratio +#' @returns A validated list of control options for log-likelihood profiling +#' @family Profiling +#' @seealso [profileLlp()] +#' @export +llpControl <- function(ofvIncrease = qchisq(0.95, df = 1), + rseTheta = 30, + itermax = 10, + ofvtol = 0.005, + paramDigits = 3, + extrapolateExpand = 1.5) { + ret <- + list( + ofvIncrease = checkmate::assert_number(ofvIncrease, lower = 0, finite = TRUE, null.ok = FALSE, na.ok = FALSE), + rseTheta = checkmate::assert_number(rseTheta, lower = 0, finite = TRUE, null.ok = FALSE, na.ok = FALSE), + ofvtol = checkmate::assert_number(ofvtol, na.ok = FALSE, lower = 1e-10, upper = 1, finite = TRUE, null.ok = FALSE), + itermax = checkmate::assert_integerish(itermax, lower = 1, any.missing = FALSE, len = 1, null.ok = FALSE, coerce = TRUE), + paramDigits = checkmate::assert_integerish(paramDigits, lower = 1, upper = 10, any.missing = FALSE, len = 1, null.ok = FALSE, coerce = TRUE), + extrapolateExpand = checkmate::assert_number(extrapolateExpand, lower = 1.001, na.ok = FALSE, finite = TRUE, null.ok = FALSE) + ) + class(ret) <- "llpControl" + ret +} + +#' @export +rxUiDeparse.llpControl <- function(object, var) { + .default <- llpControl() + .w <- nlmixr2est::.deparseDifferent(.default, object, "genRxControl") + nlmixr2est::.deparseFinal(.default, object, .w, var) +} + + +# Provide the initial estimates going up and down from the initial value +profileNlmixr2FitDataEstInitial <- function(estimates, which, ofvIncrease, rseTheta, lower, upper) { + checkmate::assert_data_frame(estimates, nrows = 1) + checkmate::assert_character(which, any.missing = FALSE, len = 1) + checkmate::assert_subset(which, choices = names(estimates)) + # use default rseTheta if not available + if (which %in% names(rseTheta)) { + currentRseTheta <- rseTheta[which] + } else { + currentRseTheta <- 30 + } + + # The -1,1 makes the estimate go up and down + ret <- estimates[[which]] + c(-1, 1)*ofvIncrease*currentRseTheta/100*abs(estimates[[which]]) + # Ensure that the estimate is within the bounds with a slight margin + margin <- sqrt(.Machine$double.eps) + pmax(pmin(ret, upper - margin), lower + margin) +} + +# Support functions ---- + +messageProfileComplete <- function(which, direction, msg) { + message( + "Profiling ", + which, + " ", + ifelse(direction > 0, "above the point estimate", "below the point estimate"), + " is ", + msg + ) +} diff --git a/R/reexport.R b/R/reexport.R index 3cc8e1c..1a3ae31 100644 --- a/R/reexport.R +++ b/R/reexport.R @@ -21,3 +21,7 @@ nlmixr2est::nlmixr2 #' @importFrom nlmixr2est setCov nlmixr2est::setCov + +#' @importFrom rxode2 rxUiDeparse +#' @export +rxode2::rxUiDeparse diff --git a/R/utils.R b/R/utils.R new file mode 100644 index 0000000..449572c --- /dev/null +++ b/R/utils.R @@ -0,0 +1,14 @@ +#' Make a control step that is quieter and faster +#' +#' @param ctl the control object +#' @return A faster and quieter control object +#' @noRd +setQuietFastControl <- function(ctl) { + # make estimation steps quieter + ctl$print <- 0L + # make estimation steps faster + ctl$covMethod <- 0L + ctl$calcTables <- FALSE + ctl$compress <- FALSE + ctl +} diff --git a/README.Rmd b/README.Rmd index aede331..fd27275 100644 --- a/README.Rmd +++ b/README.Rmd @@ -16,8 +16,8 @@ knitr::opts_chunk$set( # nlmixr2extra -[![R-CMD-check](https://github.com/nlmixr2/nlmixr2extra/workflows/R-CMD-check/badge.svg)](https://github.com/nlmixr2/nlmixr2extra/actions) -[![Codecov test coverage](https://codecov.io/gh/nlmixr2/nlmixr2extra/branch/main/graph/badge.svg)](https://app.codecov.io/gh/nlmixr2/nlmixr2extra?branch=main) +[![R-CMD-check](https://github.com/nlmixr2/nlmixr2extra/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/nlmixr2/nlmixr2extra/actions/workflows/R-CMD-check.yaml) +[![Codecov test coverage](https://codecov.io/gh/nlmixr2/nlmixr2extra/graph/badge.svg)](https://app.codecov.io/gh/nlmixr2/nlmixr2extra) [![CRAN version](http://www.r-pkg.org/badges/version/nlmixr2extra)](https://cran.r-project.org/package=nlmixr2extra) [![CRAN total downloads](https://cranlogs.r-pkg.org/badges/grand-total/nlmixr2extra)](https://cran.r-project.org/package=nlmixr2extra) [![CRAN total downloads](https://cranlogs.r-pkg.org/badges/nlmixr2extra)](https://cran.r-project.org/package=nlmixr2extra) diff --git a/man/fixedControl.Rd b/man/fixedControl.Rd new file mode 100644 index 0000000..70b50b2 --- /dev/null +++ b/man/fixedControl.Rd @@ -0,0 +1,26 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{fixedControl} +\alias{fixedControl} +\title{Control options for fixed-value likelihood profiling} +\usage{ +fixedControl() +} +\value{ +A validated list of control options for fixed-value likelihood +profiling +} +\description{ +Control options for fixed-value likelihood profiling +} +\seealso{ +\code{\link[=profileFixed]{profileFixed()}} + +Other Profiling: +\code{\link{llpControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileFixed}()}, +\code{\link{profileLlp}()}, +\code{\link{profileNlmixr2FitCoreRet}()} +} +\concept{Profiling} diff --git a/man/llpControl.Rd b/man/llpControl.Rd new file mode 100644 index 0000000..966406c --- /dev/null +++ b/man/llpControl.Rd @@ -0,0 +1,55 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{llpControl} +\alias{llpControl} +\title{Control options for log-likelihood profiling} +\usage{ +llpControl( + ofvIncrease = qchisq(0.95, df = 1), + rseTheta = 30, + itermax = 10, + ofvtol = 0.005, + paramDigits = 3, + extrapolateExpand = 1.5 +) +} +\arguments{ +\item{ofvIncrease}{The targetted change in objective function value (3.84 +corresponds to a Chi-squared test with a 95\% confidence interval)} + +\item{rseTheta}{The relative standard error (percent) for the model +parameters. It can be missing (the default) in which case a default value +of 30\% will be applied. If given as a single number, it will be applied to +all parameters. If given as a named vector of numbers, it will be applied +to each named parameter.} + +\item{itermax}{Maximum number of likelihood profiling iterations for each +bound estimated} + +\item{ofvtol}{The relative tolerance for the objective function being close +enough to the \code{ofvIncrease}.} + +\item{paramDigits}{The number of significant digits required for the +parameter. When interpolation attempts to get smaller than that number of +significant digits, it will stop.} + +\item{extrapolateExpand}{When extrapolating outside the range previously +tested, how far should the step occur as a ratio} +} +\value{ +A validated list of control options for log-likelihood profiling +} +\description{ +Control options for log-likelihood profiling +} +\seealso{ +\code{\link[=profileLlp]{profileLlp()}} + +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileFixed}()}, +\code{\link{profileLlp}()}, +\code{\link{profileNlmixr2FitCoreRet}()} +} +\concept{Profiling} diff --git a/man/profile.nlmixr2FitCore.Rd b/man/profile.nlmixr2FitCore.Rd new file mode 100644 index 0000000..d080669 --- /dev/null +++ b/man/profile.nlmixr2FitCore.Rd @@ -0,0 +1,94 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{profile.nlmixr2FitCore} +\alias{profile.nlmixr2FitCore} +\title{Perform likelihood profiling on nlmixr2 focei fits} +\usage{ +\method{profile}{nlmixr2FitCore}( + fitted, + ..., + which = NULL, + method = c("llp", "fixed"), + control = list() +) +} +\arguments{ +\item{fitted}{The fit model} + +\item{...}{ignored} + +\item{which}{The parameter names to perform likelihood profiling on +(\code{NULL} indicates all parameters)} + +\item{method}{Method to use for profiling (see the details)} + +\item{control}{Control arguments for the \code{method}} +} +\value{ +A data.frame with one column named \code{Parameter} indicating the +parameter being fixed on that row, one column for the \code{OFV} indicating the +OFV estimated for the model at that step, one column named \code{profileBound} +indicating the estimated value for the profile likelihood and its step +above the minimum profile likelihood value, and columns for each parameter +estimate (or fixed) in the model. +} +\description{ +Perform likelihood profiling on nlmixr2 focei fits +} +\section{Log-likelihood profiling}{ +\code{method = "llp"} + +The search will stop when either the OFV is within \code{ofvtol} of the desired +OFV change or when the parameter is interpolating to more significant digits +than specified in \code{paramDigits}. The "llp" method uses the \code{profileLlp()} +function. See its help for more details. +} + +\section{Fixed points}{ +\code{method = "fixed"} + +Estimate the OFV for specific fixed values. The "fixed" method uses the +\code{profileFixed()} function. See its help for more details. +} + +\examples{ +\dontrun{ +# Likelihood profiling takes a long time to run each model multiple times, so +# be aware that running this example may take a few minutes. +oneCmt <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka ~ 0.6 + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + v <- exp(tv) + cp <- linCmt() + cp ~ add(add.sd) + }) +} + +fit <- + nlmixr2( + oneCmt, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + ) +# profile all parameters +profall <- profile(fit) + +# profile a single parameter +proftka <- profile(fit, which = "tka") +} +} +\seealso{ +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{llpControl}()}, +\code{\link{profileFixed}()}, +\code{\link{profileLlp}()}, +\code{\link{profileNlmixr2FitCoreRet}()} +} +\concept{Profiling} diff --git a/man/profileFixed.Rd b/man/profileFixed.Rd new file mode 100644 index 0000000..a04c83a --- /dev/null +++ b/man/profileFixed.Rd @@ -0,0 +1,51 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{profileFixed} +\alias{profileFixed} +\alias{profileFixedSingle} +\title{Estimate the objective function values for a model while fixing defined +parameter values} +\usage{ +profileFixed(fitted, which, control = list()) + +profileFixedSingle(fitted, which) +} +\arguments{ +\item{fitted}{The fit model} + +\item{which}{A data.frame with column names of parameters to fix and values +of the fitted value to fix (one row only).} + +\item{control}{A list passed to \code{fixedControl()} (currently unused)} +} +\value{ +\code{which} with a column named \code{OFV} added with the objective function +value of the fitted estimate fixing the parameters in the other columns +} +\description{ +Estimate the objective function values for a model while fixing defined +parameter values +} +\section{Functions}{ +\itemize{ +\item \code{profileFixedSingle()}: Estimate the objective function value for a model +while fixing a single set of defined parameter values (for use in parameter +profiling) + +}} +\seealso{ +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{llpControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileLlp}()}, +\code{\link{profileNlmixr2FitCoreRet}()} + +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{llpControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileLlp}()}, +\code{\link{profileNlmixr2FitCoreRet}()} +} +\concept{Profiling} diff --git a/man/profileLlp.Rd b/man/profileLlp.Rd new file mode 100644 index 0000000..6990120 --- /dev/null +++ b/man/profileLlp.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{profileLlp} +\alias{profileLlp} +\title{Profile confidence intervals with log-likelihood profiling} +\usage{ +profileLlp(fitted, which, control) +} +\arguments{ +\item{fitted}{The fit model} + +\item{which}{Either \code{NULL} to profile all parameters or a character vector of +parameters to estimate} + +\item{control}{A list passed to \code{llpControl()}} +} +\value{ +A data.frame with columns named "Parameter" (the parameter name(s) +that were fixed), OFV (the objective function value), and the current +estimate for each of the parameters. In addition, if any boundary is +found, the OFV increase will be indicated by the absolute value of the +"profileBound" column and if that boundary is the upper or lower boundary +will be indicated by the "profileBound" column being positive or negative, +respectively. +} +\description{ +Profile confidence intervals with log-likelihood profiling +} +\seealso{ +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{llpControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileFixed}()}, +\code{\link{profileNlmixr2FitCoreRet}()} +} +\concept{Profiling} diff --git a/man/profileNlmixr2FitCoreRet.Rd b/man/profileNlmixr2FitCoreRet.Rd new file mode 100644 index 0000000..01100f4 --- /dev/null +++ b/man/profileNlmixr2FitCoreRet.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/profile.R +\name{profileNlmixr2FitCoreRet} +\alias{profileNlmixr2FitCoreRet} +\title{Give the output data.frame for a single model for profile.nlmixr2FitCore} +\usage{ +profileNlmixr2FitCoreRet(fitted, which, fixedVal) +} +\arguments{ +\item{fitted}{The fit model} + +\item{which}{The parameter names to perform likelihood profiling on +(\code{NULL} indicates all parameters)} + +\item{fixedVal}{The value that \code{which} is fixed to in case the model does not +converge.} +} +\value{ +A data.frame with columns named "Parameter" (the parameter name(s) +that were fixed), OFV (the objective function value), and the current +estimate for each of the parameters. Omega values are given as their +variances and covariances. +} +\description{ +Give the output data.frame for a single model for profile.nlmixr2FitCore +} +\seealso{ +Other Profiling: +\code{\link{fixedControl}()}, +\code{\link{llpControl}()}, +\code{\link{profile.nlmixr2FitCore}()}, +\code{\link{profileFixed}()}, +\code{\link{profileLlp}()} +} +\concept{Profiling} diff --git a/man/reexports.Rd b/man/reexports.Rd index 78e68d7..1cc2959 100644 --- a/man/reexports.Rd +++ b/man/reexports.Rd @@ -10,6 +10,7 @@ \alias{nlmixr} \alias{nlmixr2} \alias{setCov} +\alias{rxUiDeparse} \title{Objects exported from other packages} \keyword{internal} \description{ @@ -21,6 +22,6 @@ below to see their documentation. \item{nlmixr2est}{\code{\link[nlmixr2est:nlmixr2]{nlmixr}}, \code{\link[nlmixr2est]{nlmixr2}}, \code{\link[nlmixr2est]{nlmixrWithTiming}}, \code{\link[nlmixr2est]{setCov}}} - \item{rxode2}{\code{\link[rxode2]{ini}}, \code{\link[rxode2]{model}}} + \item{rxode2}{\code{\link[rxode2]{ini}}, \code{\link[rxode2]{model}}, \code{\link[rxode2]{rxUiDeparse}}} }} diff --git a/tests/testthat/test-profile.R b/tests/testthat/test-profile.R new file mode 100644 index 0000000..27126ae --- /dev/null +++ b/tests/testthat/test-profile.R @@ -0,0 +1,230 @@ +test_that("profileNlmixr2FitDataEstInitial", { + # Must have one-row input + expect_error( + profileNlmixr2FitDataEstInitial(estimates = data.frame(A = 1:2)) + ) + + expect_equal( + profileNlmixr2FitDataEstInitial( + estimates = data.frame(A = 1), + which = "A", + ofvIncrease = 1.92, + rseTheta = c(A=100), + lower = -100, upper = 200 + ), + c(-0.92, 2.92) + ) + # Bounds are respected + expect_equal( + profileNlmixr2FitDataEstInitial( + estimates = data.frame(A = 1), + which = "A", + ofvIncrease = 1.92, + rseTheta = c(A=100), + lower = 0, upper = 200 + ), + c(sqrt(.Machine$double.eps), 2.92) + ) +}) + +test_that("profileNlmixr2FitCoreRet", { + # Variance and covariance is correctly captured + one.compartment <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka ~ 0.6 + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + v <- exp(tv) + cp <- linCmt() + cp ~ add(add.sd) + }) + } + + fit <- + suppressMessages(nlmixr2( + one.compartment, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + )) + withoutCov <- profileNlmixr2FitCoreRet(fit, which = "tka") + expect_s3_class(withoutCov, "data.frame") + expect_named(withoutCov, expected = c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka")) + + one.compartment <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka + eta.cl ~ c(0.6, 0.1, 0.2) + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl + eta.cl) + v <- exp(tv) + cp <- linCmt() + cp ~ add(add.sd) + }) + } + + fit <- + suppressMessages(nlmixr2( + one.compartment, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + )) + withCov <- profileNlmixr2FitCoreRet(fit, which = "tka") + expect_s3_class(withCov, "data.frame") + expect_named(withCov, expected = c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka", "eta.cl", "cov(eta.cl,eta.ka)")) +}) + +test_that("profileFixed", { + # fix most of the parameters so that it estimates faster + one.compartment <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka ~ 0.6 + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + v <- exp(tv) + cp <- linCmt() + cp ~ add(add.sd) + }) + } + + fit <- + suppressMessages(nlmixr2( + one.compartment, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + )) + + testFixed <- + suppressMessages( + profile(fit, which = data.frame(tka = log(c(1.4, 1.6, 1.8))), method = "fixed") + ) + expect_s3_class(testFixed, "data.frame") + expect_named(testFixed, expected = c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka")) + expect_equal(nrow(testFixed), 3) + + # Fix multiple parameters simultaneously + testFixedMulti <- + suppressMessages( + profile( + fit, + which = + data.frame( + tka = log(c(1.4, 1.6, 1.8)), + tcl = log(c(2.6, 2.7, 2.8)) + ), + method = "fixed" + ) + ) + expect_s3_class(testFixedMulti, "data.frame") + expect_named(testFixedMulti, expected = c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka")) + expect_equal(nrow(testFixedMulti), 3) + expect_equal(testFixedMulti$Parameter, rep("tka,tcl", 3)) +}) + +test_that("profile a standard model", { + # fix most of the parameters so that it estimates faster + one.compartment <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka ~ 0.6 + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl) + v <- exp(tv) + cp <- linCmt() + cp ~ add(add.sd) + }) + } + + fit <- + suppressMessages(nlmixr2( + one.compartment, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + )) + + # All parameters + profall <- suppressMessages(profile(fit)) + expect_s3_class(profall, "data.frame") + expect_named(profall, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "profileBound")) + + # A single parameter + proftka <- suppressMessages(profile(fit, which = "tka")) + expect_s3_class(proftka, "data.frame") + expect_named(proftka, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "profileBound")) + + # A fixed parameter + expect_warning( + proftv <- profile(fit, which = "tv"), + regexp = "OFV decreased while profiling" + ) + expect_s3_class(proftv, "data.frame") + expect_named(proftv, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "profileBound")) + + # Residual error + profadd.sd <- profile(fit, which = "add.sd") + expect_s3_class(profadd.sd, "data.frame") + expect_named(profadd.sd, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "profileBound")) +}) + +test_that("profile a standard model with correlated etas", { + # fix most of the parameters so that it estimates faster + one.compartment <- function() { + ini({ + tka <- log(1.57) + tcl <- log(2.72) + tv <- fixed(log(31.5)) + eta.ka ~ 0.6 + eta.cl ~ 0.1 + eta.v ~ 0.2 + add.sd <- 0.7 + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl + eta.cl) + v <- exp(tv + eta.v) + cp <- linCmt() + cp ~ add(add.sd) + }) + } + + fit <- + suppressMessages(nlmixr2( + one.compartment, data = nlmixr2data::theo_sd, est="focei", control = list(print=0) + )) + + # All parameters + profall <- suppressMessages(profile(fit)) + expect_s3_class(profall, "data.frame") + expect_named(profall, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka", "profileBound")) + + # A single parameter + proftka <- suppressMessages(profile(fit, which = "tka")) + expect_s3_class(proftka, "data.frame") + expect_named(proftka, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka", "profileBound")) + + # A fixed parameter + expect_warning( + proftv <- profile(fit, which = "tv"), + regexp = "OFV decreased while profiling" + ) + expect_s3_class(proftv, "data.frame") + expect_named(proftv, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka", "profileBound")) + + # Residual error + profadd.sd <- profile(fit, which = "add.sd") + expect_s3_class(profadd.sd, "data.frame") + expect_named(profadd.sd, c("Parameter", "OFV", "tka", "tcl", "tv", "add.sd", "eta.ka", "profileBound")) +})