From 4ab2688fd237c3c7266a85c15a14128712ec1920 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 15:46:42 -0500 Subject: [PATCH 01/11] Add .rmDiag(ui, diag) in preparation for %>% ini(diag) --- R/piping-ini.R | 103 ++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 101 insertions(+), 2 deletions(-) diff --git a/R/piping-ini.R b/R/piping-ini.R index b9b8a0744..8d7d1836e 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -1,3 +1,11 @@ +#' Message about fixing or unfixing a parameter +#' +#' @param ini this is the iniDf data frame +#' @param w this indicates the row number of the item that is fixed or +#' unfixed +#' @param fixedValue this is a boolean +#' @noRd +#' @author Matthew L. Fidler .msgFix<- function(ini, w, fixedValue) { lapply(w, function(.w) { if (ini$fix[.w] != fixedValue) { @@ -10,6 +18,18 @@ }) } +#' This modifies the iniDf to fix (or unfix) parameters and related +#' values +#' +#' Note that the block of etas will be fixed/unfixed when a single +#' value is fixed/unfixed +#' +#' @param ini iniDf data.frame +#' @param w which item will be fixed +#' @param fixedValue should this be fixed `TRUE` or unfixed `FALSE` +#' @return nothing, called for side effects +#' @noRd +#' @author Matthew L. Fidler .iniModifyFixedForThetaOrEtablock <- function(ini, w, fixedValue) { if (rxode2.verbose.pipe) { .msgFix(ini, w, fixedValue) @@ -215,9 +235,9 @@ rbind(ini,.ini2) } -#' This function handles the lotri process and integrates into current UI +#' This function handles the lotri process and integrates into current UI #' -#' This will update the matrix and integrate the initial estimates in the UI +#' This will update the matrix and integrate the initial estimates in the UI #' #' @param mat Lotri processed matrix from the piping ini function #' @@ -892,3 +912,82 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { ini(.ret) <- iniDf .ret } + +#' This removes the off-diagonal BSV from a rxode2 model +#' +#' @param ui rxode2 ui model +#' +#' @param diag character vector of diagonal values to remove +#' +#' @return model with off-diagonals removed +#' +#' @export +#' +#' @author Matthew L. Fidler +#' +#' @examples +#' +#' one.compartment <- function() { +#' ini({ +#' tka <- log(1.57); label("Ka") +#' tcl <- log(2.72); label("Cl") +#' tv <- log(31.5); label("V") +#' eta.ka ~ c(0.6) +#' eta.cl ~ c(0.01, 0.3) +#' eta.v ~ c(0.01, 0.01, 0.1) +#' add.sd <- 0.7 +#' }) +#' model({ +#' ka <- exp(tka + eta.ka) +#' cl <- exp(tcl + eta.cl) +#' v <- exp(tv + eta.v) +#' d/dt(depot) = -ka * depot +#' d/dt(center) = ka * depot - cl / v * center +#' cp = center / v +#' cp ~ add(add.sd) +#' }) +#' } +#' +.rmDiag <- function(ui, diag=character(0)) { + .ui <- rxode2::assertRxUi(ui) + .iniDf <- .ui$iniDf + .theta <- .iniDf[!is.na(.iniDf$ntheta),,drop=FALSE] + .eta <- .iniDf[is.na(.iniDf$ntheta),,drop=FALSE] + if (length(diag) == 0) { + .w <- which(.eta$neta1 == .eta$neta2) + .rmNames <- .eta[-.w, "name"] + .eta <- .eta[which(.w),, drop=FALSE] + .iniDf <- rbind(.theta, .eta) + ini(.ui) <- .iniDf + } else { + .rmNames <- character(0) + for (.e in diag) { + .w <- which(.eta$name == .e) + if (length(.w) == 1L) { + .n <- .eta$neta1[.w] + .w <- vapply(seq_along(.eta$neta1), + function(i) { + if (.eta$neta1[i] == .eta$neta2[i]) { + TRUE + } else if (.eta$neta1[i] == .n && .eta$neta2[i] != .n) { + FALSE + } else if (.eta$neta2[i] == .n && .eta$neta1[i] != .n) { + FALSE + } else { + TRUE + } + }, logical(1), USE.NAMES = TRUE) + .rmNames <- c(.rmNames, .eta$name[!.w]) + .eta <- .eta[.w,,drop=FALSE] + } + } + .iniDf <- rbind(.theta, .eta) + ini(.ui) <- .iniDf + } + if (rxode2.verbose.pipe) { + for (.v in .rmNames) { + .minfo(paste0("remove covariance {.code ", .v, "}")) + } + } + .ui +} From ab146a1226964d251b18ce9d31d10703adcddf8e Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 17:33:09 -0500 Subject: [PATCH 02/11] Add diag() handling --- R/piping-ini.R | 106 +++++++++++++++++++++++++++++++------------------ R/piping.R | 8 ++++ R/ui.R | 25 ++++++++---- 3 files changed, 93 insertions(+), 46 deletions(-) diff --git a/R/piping-ini.R b/R/piping-ini.R index 8d7d1836e..fb1213469 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -448,7 +448,6 @@ # This likely cannot be reached because all scenarios should be handled # above in the input checking. The line remains in the code defensively. stop("Cannot find parameter '", append, "'", call.=FALSE) # nocov - } else if (appendClean == wLhs) { warning("parameter '", lhs, "' set to be moved after itself, no change in order made", call. = FALSE) @@ -608,6 +607,17 @@ #' @keywords internal #' @export .iniHandleLine <- function(expr, rxui, envir=parent.frame(), append = NULL) { + if (.matchesLangTemplate(expr, str2lang("~diag()"))) { + .iniHandleDiag(expr=NULL, rxui=rxui) + return(invisible()) + } else if (length(expr) == 2L && + identical(expr[[1]], quote(`~`)) && + is.call(expr[[2]]) && length(expr[[2]]) >= 2L && + identical(expr[[2]][[1]], quote(`diag`))) { + # .matchesLangTemplate(expr, str2lang("~diag(.)")) doesn't work + .iniHandleDiag(expr=expr, rxui=rxui) + return(invisible()) + } # Convert all variations on fix, fixed, FIX, FIXED; unfix, unfixed, UNFIX, # UNFIXED to fix and unfix to simplify all downstream operations expr <- .iniSimplifyFixUnfix(expr) @@ -627,7 +637,6 @@ } else if (.matchesLangTemplate(expr, str2lang("unfix(.name)"))) { expr <- as.call(list(quote(`<-`), expr[[2]], quote(`unfix`))) } - if (.matchesLangTemplate(expr, str2lang(".name <- label(.)"))) { .iniHandleLabel(expr=expr, rxui=rxui, envir=envir) } else if (.matchesLangTemplate(expr, str2lang(".name <- backTransform(.)"))) { @@ -913,52 +922,24 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { .ret } -#' This removes the off-diagonal BSV from a rxode2 model +#' This removes the off-diagonal BSV from a rxode2 iniDf #' #' @param ui rxode2 ui model #' #' @param diag character vector of diagonal values to remove #' -#' @return model with off-diagonals removed -#' -#' @export -#' +#' @return iniDf with modified diagonal +#' @noRd #' @author Matthew L. Fidler -#' -#' @examples -#' -#' one.compartment <- function() { -#' ini({ -#' tka <- log(1.57); label("Ka") -#' tcl <- log(2.72); label("Cl") -#' tv <- log(31.5); label("V") -#' eta.ka ~ c(0.6) -#' eta.cl ~ c(0.01, 0.3) -#' eta.v ~ c(0.01, 0.01, 0.1) -#' add.sd <- 0.7 -#' }) -#' model({ -#' ka <- exp(tka + eta.ka) -#' cl <- exp(tcl + eta.cl) -#' v <- exp(tv + eta.v) -#' d/dt(depot) = -ka * depot -#' d/dt(center) = ka * depot - cl / v * center -#' cp = center / v -#' cp ~ add(add.sd) -#' }) -#' } -#' -.rmDiag <- function(ui, diag=character(0)) { - .ui <- rxode2::assertRxUi(ui) - .iniDf <- .ui$iniDf +.iniDfRmDiag <- function(iniDf, diag=character(0)) { + .iniDf <- iniDf .theta <- .iniDf[!is.na(.iniDf$ntheta),,drop=FALSE] .eta <- .iniDf[is.na(.iniDf$ntheta),,drop=FALSE] if (length(diag) == 0) { .w <- which(.eta$neta1 == .eta$neta2) .rmNames <- .eta[-.w, "name"] - .eta <- .eta[which(.w),, drop=FALSE] + .eta <- .eta[.w,, drop=FALSE] .iniDf <- rbind(.theta, .eta) - ini(.ui) <- .iniDf } else { .rmNames <- character(0) for (.e in diag) { @@ -981,13 +962,62 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { .eta <- .eta[.w,,drop=FALSE] } } + .mat <- lotri::as.lotri(.eta) + .mat <- lotri::rcm(.mat) + class(.mat) <- c("lotriFix", class(.mat)) + .eta <- as.data.frame(.mat) + .eta$err <- NA_character_ .iniDf <- rbind(.theta, .eta) - ini(.ui) <- .iniDf } if (rxode2.verbose.pipe) { for (.v in .rmNames) { .minfo(paste0("remove covariance {.code ", .v, "}")) } } - .ui + .iniDf +} + +.iniHandleDiag <- function(expr, rxui){ + if (is.null(expr)) { + assign("iniDf", .iniDfRmDiag(rxui$iniDf), envir=rxui) + } else { + # now get the variables in the diag expression + .env <- new.env(parent=emptyenv()) + .env$names <- character(0) + .f <- function(x) { + if (is.name(x)) { + .env$names <- c(.env$names, as.character(x)) + } else if (is.call(x)) { + lapply(lapply(seq_along(x)[-1], function(i) {x[[i]]}), .f) + } + } + expr <- expr[[2]] + lapply(seq_along(expr)[-1], + function(i) { + .f(expr[[i]]) + }) + ## .mat <- as.matrix(lotri::as.lotri(rxui$iniDf)) + ## attr(.mat, "lotriEst") <- NULL + ## class(.mat) <- NULL + ## .tmp <- expand.grid(r=.env$names, c=.env$names) + ## .tmp <- .tmp[.tmp$r != .tmp$c,] + ## .tmp$r <- paste(.tmp$r) + ## .tmp$c <- paste(.tmp$c) + ## .tmp$r <- vapply(.tmp$r, + ## function(x) { + ## which(x == colnames(.mat)) + ## }, integer(1)) + ## .tmp$c <- vapply(.tmp$c, + ## function(x) { + ## which(x == colnames(.mat)) + ## }, integer(1)) + ## for (.i in seq_along(.tmp$r)) { + ## .mat[.tmp$r[.i], .tmp$c[.i]] <- 0.0 + ## } + ## .mat <- lotri::rcm(.mat) + ## class(.mat) <- c("lotriFix", class(.mat)) + ## .ini <- as.data.frame(.mat) + # remove covariates between the expressions + assign("iniDf", .iniDfRmDiag(rxui$iniDf, .env$names), envir=rxui) + } } diff --git a/R/piping.R b/R/piping.R index 099c1ee69..e59bcadd9 100644 --- a/R/piping.R +++ b/R/piping.R @@ -373,9 +373,17 @@ # Capture empty arguments (rxode2#688) warning("empty argument ignored") return(NULL) + } else if (is.symbol(.quoted) && + identical(.quoted, quote(`diag`))) { + .quoted <- str2lang("~diag()") + } else if (length(.quoted) >= 1 && + identical(.quoted[[1]], quote(`diag`))) { + .quoted <- as.call(c(list(quote(`~`)), .quoted)) } else if (length(.quoted) == 1) { .bracket[i] <- TRUE assign(".bracket", .bracket, envir=.env) + } else if (identical(.quoted[[1]], quote(`diag`))) { + } else if (identical(.quoted[[1]], quote(`{`)) || identical(.quoted[[1]], quote(`c`)) || identical(.quoted[[1]], quote(`list`))) { diff --git a/R/ui.R b/R/ui.R index da481d9cb..d77738e47 100644 --- a/R/ui.R +++ b/R/ui.R @@ -137,14 +137,17 @@ #' #' 'omega' values can be set as a single value or as the values of a #' lower-triangular matrix. The values may be set as either a -#' variance-covariance matrix (the default) or as a correlation matrix for the -#' off-diagonals with the standard deviations on the diagonals. Names may be -#' set on the left side of the \code{~}. To set a variance-covariance matrix -#' with variance values of 2 and 3 and a covariance of -2.5 use \code{~c(2, 2.5, -#' 3)}. To set the same matrix with names of \code{iivKa} and \code{iivCL}, use -#' \code{iivKa + iivCL~c(2, 2.5, 3)}. To set a correlation matrix with standard -#' deviations on the diagonal, use \code{cor()} like \code{iivKa + iivCL~cor(2, -#' -0.5, 3)}. +#' variance-covariance matrix (the default) or as a correlation matrix +#' for the off-diagonals with the standard deviations on the +#' diagonals. Names may be set on the left side of the \code{~}. To +#' set a variance-covariance matrix with variance values of 2 and 3 +#' and a covariance of -2.5 use \code{~c(2, 2.5, 3)}. To set the same +#' matrix with names of \code{iivKa} and \code{iivCL}, use \code{iivKa +#' + iivCL~c(2, 2.5, 3)}. To set a correlation matrix with standard +#' deviations on the diagonal, use \code{cor()} like \code{iivKa + +#' iivCL~cor(2, -0.5, 3)}. As of rxode2 3.0 you can also use +#' \code{iivKa ~ 2, iivCL ~ c(2.5, 3)} for covariance matrices as +#' well. #' #' Values may be fixed (and therefore not estimated) using either the name #' \code{fixed} at the end of the assignment or by calling \code{fixed()} as a @@ -173,6 +176,12 @@ #' has a label of "Typical Value of Clearance (L/hr)" is \code{tvCL <- 1; #' label("Typical Value of Clearance (L/hr)")}. #' +#' Off diagonal values of 'omega' can be set to zero using the +#' \code{diag()} or \code{diag(iivKa, iivCL)} for example removing all +#' off-diagonals can be removed with `ini(diag())`, or the off +#' diagonals between clearance and ka could be removed by +#' \code{ini(diag(iivKa, iivCL)}. +#' #' \code{rxode2}/\code{nlmixr2} will attempt to determine some #' back-transformations for the user. For example, \code{CL <- exp(tvCL)} will #' detect that \code{tvCL} must be back-transformed by \code{exp()} for easier From b60d77994f11e480c2a21c690bcbe11a0d8bdb21 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 17:50:30 -0500 Subject: [PATCH 03/11] Add some tests for diag handling --- tests/testthat/test-piping-ini.R | 69 ++++++++++++++++++++++++++++++++ 1 file changed, 69 insertions(+) diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R index 38038204f..558d70274 100644 --- a/tests/testthat/test-piping-ini.R +++ b/tests/testthat/test-piping-ini.R @@ -908,3 +908,72 @@ test_that("parameters can be promoted from covariate to parameter with bounds (# fixed = TRUE ) }) + +test_that("ini(diag) tests", { + + mod2 <- function() { + ini({ + lka ~ 0.45 + lcl ~ c(0.01, 1) + lvc ~ c(-0.01, 0.01, 3.45) + lfun ~ c(-0.1, 0.1, 0.01, 4) + }) + model({ + ka <- exp(lka) + cl <- exp(lcl) + vc <- exp(lvc) + kel <- cl / vc + d/dt(depot) <- -ka*depot + d/dt(central) <- ka*depot-kel*central + cp <- central / vc + lfun + }) + } + + # Will reorder + tmp <- mod2 %>% ini(diag(lcl, lvc)) + + expect_equal(tmp$omega, + lotri({ + lfun ~ 4 + lka ~ c(-0.1, 0.45) + lvc ~ 3.45 + lcl ~ 1 + })) + + tmp <- mod2 %>% ini(diag) + expect_equal(tmp$omega, + lotri({ + lka ~ 0.45 + lcl ~ 1 + lvc ~ 3.45 + lfun ~ 4 + })) + + mod <- function() { + ini({ + lka ~ 0.45 + lcl ~ c(0.01, 1) + lvc ~ c(-0.01, 0.01, 3.45) + }) + model({ + ka <- exp(lka) + cl <- exp(lcl) + vc <- exp(lvc) + kel <- cl / vc + d/dt(depot) <- -ka*depot + d/dt(central) <- ka*depot-kel*central + cp <- central / vc + }) + } + + + tmp <- mod %>% ini(diag) + + expect_equal(tmp$omega, + lotri({ + lka ~ 0.45 + lcl ~ 1 + lvc ~ 3.45 + })) + +}) From 3b29dd5f8c5103b3c017b6dc55177a076143ecc0 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 17:54:56 -0500 Subject: [PATCH 04/11] Add to news --- NEWS.md | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/NEWS.md b/NEWS.md index d5496bbfb..4f07d6774 100644 --- a/NEWS.md +++ b/NEWS.md @@ -49,6 +49,12 @@ ## New features +- You can remove covariances for every omega by piping with `%>% + ini(diag())` you can be a bit more granular by removing all + covariances that have either `eta.ka` or `eta.cl` by: `%>% + ini(diag(eta.ka, eta.cl))` or anything with correlations with + `eta.cl` with `%>% ini(diag(eta.cl))` + - You can specify the type of interpolation applied for added dosing records (or other added records) for columns that are kept with the `keep=` option in `rxSolve()`. This new option is From fa9a742bf257a18a39632079686e9c5f0197a748 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 17:57:35 -0500 Subject: [PATCH 05/11] Add single item test --- tests/testthat/test-piping-ini.R | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R index 558d70274..69738b3e1 100644 --- a/tests/testthat/test-piping-ini.R +++ b/tests/testthat/test-piping-ini.R @@ -949,6 +949,16 @@ test_that("ini(diag) tests", { lfun ~ 4 })) + tmp <- mod2 %>% ini(diag(lvc)) + + expect_equal(tmp$omega, + lotri({ + lfun ~ 4 + lcl ~ c(0.1, 1) + lka ~ c(-0.1, 0.01, 0.45) + lvc ~ 3.45 + })) + mod <- function() { ini({ lka ~ 0.45 From 3e491c2bc0b3c87d0c14248a03400ed892eadd51 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 18:14:50 -0500 Subject: [PATCH 06/11] Add -cov() and -cor() --- R/piping-ini.R | 45 ++++++++++++++++---------------- tests/testthat/test-piping-ini.R | 39 +++++++++++++++++++++++++-- 2 files changed, 60 insertions(+), 24 deletions(-) diff --git a/R/piping-ini.R b/R/piping-ini.R index fb1213469..5a677beb6 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -631,6 +631,13 @@ expr <- as.call(list(quote(`-`), expr[[2]])) } + # now handle dropping covariances + if (.matchesLangTemplate(expr, str2lang("-cov(.name, .name)")) || + .matchesLangTemplate(expr, str2lang("-cor(.name, .name)"))) { + .iniHandleRmCov(expr=expr, rxui=rxui) + return(invisible()) + } + # Convert fix(name) or unfix(name) to name <- fix or name <- unfix if (.matchesLangTemplate(expr, str2lang("fix(.name)"))) { expr <- as.call(list(quote(`<-`), expr[[2]], quote(`fix`))) @@ -977,6 +984,22 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { .iniDf } +.iniHandleRmCov <- function(expr, rxui) { + .iniDf <- rxui$iniDf + .theta <- .iniDf[!is.na(.iniDf$ntheta),, drop = FALSE] + .eta <- .iniDf[is.na(.iniDf$ntheta),, drop = FALSE] + .mat <- lotri::as.lotri(.eta) + .v1 <- which(as.character(expr[[2]][[2]])==dimnames(.mat)[[1]]) + .v2 <- which(as.character(expr[[2]][[3]])==dimnames(.mat)[[1]]) + .mat[.v1, .v2] <- .mat[.v2, .v1] <- 0 + .mat <- lotri::rcm(.mat) + class(.mat) <- c("lotriFix", class(.mat)) + .eta <- as.data.frame(.mat) + .eta$err <- NA_character_ + .iniDf <- rbind(.theta, .eta) + assign("iniDf", .iniDf, envir=rxui) +} + .iniHandleDiag <- function(expr, rxui){ if (is.null(expr)) { assign("iniDf", .iniDfRmDiag(rxui$iniDf), envir=rxui) @@ -996,28 +1019,6 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { function(i) { .f(expr[[i]]) }) - ## .mat <- as.matrix(lotri::as.lotri(rxui$iniDf)) - ## attr(.mat, "lotriEst") <- NULL - ## class(.mat) <- NULL - ## .tmp <- expand.grid(r=.env$names, c=.env$names) - ## .tmp <- .tmp[.tmp$r != .tmp$c,] - ## .tmp$r <- paste(.tmp$r) - ## .tmp$c <- paste(.tmp$c) - ## .tmp$r <- vapply(.tmp$r, - ## function(x) { - ## which(x == colnames(.mat)) - ## }, integer(1)) - ## .tmp$c <- vapply(.tmp$c, - ## function(x) { - ## which(x == colnames(.mat)) - ## }, integer(1)) - ## for (.i in seq_along(.tmp$r)) { - ## .mat[.tmp$r[.i], .tmp$c[.i]] <- 0.0 - ## } - ## .mat <- lotri::rcm(.mat) - ## class(.mat) <- c("lotriFix", class(.mat)) - ## .ini <- as.data.frame(.mat) - # remove covariates between the expressions assign("iniDf", .iniDfRmDiag(rxui$iniDf, .env$names), envir=rxui) } } diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R index 69738b3e1..65c7e610c 100644 --- a/tests/testthat/test-piping-ini.R +++ b/tests/testthat/test-piping-ini.R @@ -909,7 +909,7 @@ test_that("parameters can be promoted from covariate to parameter with bounds (# ) }) -test_that("ini(diag) tests", { +test_that("ini(diag) and ini(-cov()) tests", { mod2 <- function() { ini({ @@ -929,9 +929,44 @@ test_that("ini(diag) tests", { }) } + tmp <- mod2 %>% ini(-cov(lcl, lvc)) + expect_equal(tmp$omega, + lotri({ + lvc ~ 3.45 + lfun ~ c(0.01, 4) + lka ~ c(-0.01, -0.1, 0.45) + lcl ~ c(0, 0.1, 0.01, 1) + })) + + tmp <- mod2 %>% ini(-cor(lcl, lvc)) + expect_equal(tmp$omega, + lotri({ + lvc ~ 3.45 + lfun ~ c(0.01, 4) + lka ~ c(-0.01, -0.1, 0.45) + lcl ~ c(0, 0.1, 0.01, 1) + })) + + tmp <- mod2 %>% ini(cor(lcl, lvc) <- NULL) + expect_equal(tmp$omega, + lotri({ + lvc ~ 3.45 + lfun ~ c(0.01, 4) + lka ~ c(-0.01, -0.1, 0.45) + lcl ~ c(0, 0.1, 0.01, 1) + })) + + tmp <- mod2 %>% ini(cor(lcl, lvc) ~ NULL) + expect_equal(tmp$omega, + lotri({ + lvc ~ 3.45 + lfun ~ c(0.01, 4) + lka ~ c(-0.01, -0.1, 0.45) + lcl ~ c(0, 0.1, 0.01, 1) + })) + # Will reorder tmp <- mod2 %>% ini(diag(lcl, lvc)) - expect_equal(tmp$omega, lotri({ lfun ~ 4 From bd88b7577db07f1ba722310f144aa263a5e21bd0 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 18:24:04 -0500 Subject: [PATCH 07/11] Fix cov <- NULL and cov ~ NULL --- R/piping-ini.R | 9 ++++++--- R/piping.R | 2 +- tests/testthat/test-piping-ini.R | 1 + 3 files changed, 8 insertions(+), 4 deletions(-) diff --git a/R/piping-ini.R b/R/piping-ini.R index 5a677beb6..643f6f131 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -625,9 +625,12 @@ # downstream operations expr <- .iniSimplifyAssignArrow(expr) - if (.matchesLangTemplate(expr, str2lang(".name <- NULL"))) { - expr <- as.call(list(quote(`-`), expr[[2]])) - } else if (.matchesLangTemplate(expr, str2lang(".name ~ NULL"))) { + if (.matchesLangTemplate(expr, str2lang(".name <- NULL")) || + .matchesLangTemplate(expr, str2lang(".name ~ NULL")) || + .matchesLangTemplate(expr, str2lang("cov(.name, .name) <- NULL")) || + .matchesLangTemplate(expr, str2lang("cor(.name, .name) <- NULL")) || + .matchesLangTemplate(expr, str2lang("cov(.name, .name) ~ NULL")) || + .matchesLangTemplate(expr, str2lang("cor(.name, .name) ~ NULL"))) { expr <- as.call(list(quote(`-`), expr[[2]])) } diff --git a/R/piping.R b/R/piping.R index e59bcadd9..8c3beaddd 100644 --- a/R/piping.R +++ b/R/piping.R @@ -392,7 +392,7 @@ } else if (identical(.quoted[[1]], quote(`as.formula`))) { .quoted <- .quoted[[2]] } else if (identical(.quoted[[1]], quote(`~`))) { - if (length(.quoted) == 3L) { + if (length(.quoted) == 3L && !is.null(.quoted[[3]])) { .quoted[[3]] <- .iniSimplifyFixUnfix(.quoted[[3]]) if (identical(.quoted[[3]], quote(`fix`)) || identical(.quoted[[3]], quote(`unfix`))) { diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R index 65c7e610c..f330ab591 100644 --- a/tests/testthat/test-piping-ini.R +++ b/tests/testthat/test-piping-ini.R @@ -948,6 +948,7 @@ test_that("ini(diag) and ini(-cov()) tests", { })) tmp <- mod2 %>% ini(cor(lcl, lvc) <- NULL) + expect_equal(tmp$omega, lotri({ lvc ~ 3.45 From 1e3a3e260f2f9c0891e8b2e3d661dea282ead4ca Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 18:36:32 -0500 Subject: [PATCH 08/11] Add to vignette and news --- NEWS.md | 3 ++ vignettes/articles/Modifying-Models.Rmd | 45 +++++++++++++++++++++++++ 2 files changed, 48 insertions(+) diff --git a/NEWS.md b/NEWS.md index 4f07d6774..a6d351caf 100644 --- a/NEWS.md +++ b/NEWS.md @@ -55,6 +55,9 @@ ini(diag(eta.ka, eta.cl))` or anything with correlations with `eta.cl` with `%>% ini(diag(eta.cl))` +- You can also remove individual covariances by `%>% ini(-cov(a, b))` + or `%>% ini(-cor(a,b))`. + - You can specify the type of interpolation applied for added dosing records (or other added records) for columns that are kept with the `keep=` option in `rxSolve()`. This new option is diff --git a/vignettes/articles/Modifying-Models.Rmd b/vignettes/articles/Modifying-Models.Rmd index ff8fd7f00..b421ee44e 100644 --- a/vignettes/articles/Modifying-Models.Rmd +++ b/vignettes/articles/Modifying-Models.Rmd @@ -502,6 +502,8 @@ The common items you want to do with initial estimates are: - Reorder parameters +- Remove covariances between all parameters or a group of parameters + You may wish to create your own functions; we will discuss this too. @@ -830,6 +832,49 @@ f8 <- f |> print(f8) ``` +### Removing covariances between between subject varaibilities + +There are two approaches to removing covarinaces for between subject +variabilities: `diag()` and `-cov(var1, var2)` + +The `diag()` removes either all covariance elements (with no +arguments) or any covariance elements included in the argument list: + +```{r diagExamples} +fd <- function() { + ini({ + tka <- 0.45 + tcl <- 1 + tv <- 3.45 + add.sd <- c(0, 0.7) + eta.ka ~ 0.6 + eta.v ~ c(0.01, 0.1) + eta.cl ~ c(0.01, 0.01, 1) + }) + model({ + ka <- exp(tka + eta.ka) + cl <- exp(tcl + eta.cl) + v <- exp(tv + eta.v) + d/dt(depot) <- -ka * depot + d/dt(center) <- ka * depot - cl/v * center + cp <- center/v + cp ~ add(add.sd) + }) +} + +# If you want to remove all covariances you can use diag() with no +# arguments + +fd %>% ini(diag()) + +# if you want to remove only covariances with eta.ka you can use: +fd %>% ini(diag(eta.ka)) + +# if you want to remove only the covariances with eta.ka and eta.v you can use: +fd %>% ini(-cov(eta.ka, eta.v)) + +``` + ### More granular access of initial conditions Just like with `model()` you can modify the underlying data frame that From 7d537e393c4a5033ae6cbd6b2d42c8d45584f6b6 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 22:00:59 -0500 Subject: [PATCH 09/11] Add errors --- R/piping-ini.R | 18 ++++++++++++++++-- tests/testthat/test-piping-ini.R | 13 +++++++++++++ 2 files changed, 29 insertions(+), 2 deletions(-) diff --git a/R/piping-ini.R b/R/piping-ini.R index 643f6f131..226316a74 100644 --- a/R/piping-ini.R +++ b/R/piping-ini.R @@ -970,6 +970,8 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { }, logical(1), USE.NAMES = TRUE) .rmNames <- c(.rmNames, .eta$name[!.w]) .eta <- .eta[.w,,drop=FALSE] + } else { + stop("cannot find parameter '", .e, "' for covariance removal", call.=FALSE) } } .mat <- lotri::as.lotri(.eta) @@ -992,8 +994,20 @@ zeroRe <- function(object, which = c("omega", "sigma"), fix = TRUE) { .theta <- .iniDf[!is.na(.iniDf$ntheta),, drop = FALSE] .eta <- .iniDf[is.na(.iniDf$ntheta),, drop = FALSE] .mat <- lotri::as.lotri(.eta) - .v1 <- which(as.character(expr[[2]][[2]])==dimnames(.mat)[[1]]) - .v2 <- which(as.character(expr[[2]][[3]])==dimnames(.mat)[[1]]) + .n1 <- as.character(expr[[2]][[2]]) + .v1 <- which(.n1==dimnames(.mat)[[1]]) + if (length(.v1) != 1) { + stop("cannot find parameter '", .n1, "' for covariance removal", call.=FALSE) + } + .n2 <- as.character(expr[[2]][[3]]) + .v2 <- which(.n2==dimnames(.mat)[[1]]) + if (length(.v2) != 1) { + stop("cannot find parameter '", .n2, "' for covariance removal", call.=FALSE) + } + if (rxode2.verbose.pipe) { + .minfo(paste0("remove covariance {.code (", .n1, ", ", .n2, ")}")) + } + .mat[.v1, .v2] <- .mat[.v2, .v1] <- 0 .mat <- lotri::rcm(.mat) class(.mat) <- c("lotriFix", class(.mat)) diff --git a/tests/testthat/test-piping-ini.R b/tests/testthat/test-piping-ini.R index f330ab591..b327cf0fe 100644 --- a/tests/testthat/test-piping-ini.R +++ b/tests/testthat/test-piping-ini.R @@ -929,6 +929,16 @@ test_that("ini(diag) and ini(-cov()) tests", { }) } + expect_error( + mod2 %>% ini(diag(lcl, matt)), + "matt" + ) + + expect_error( + mod2 %>% ini(diag(matt, lcl)), + "matt" + ) + tmp <- mod2 %>% ini(-cov(lcl, lvc)) expect_equal(tmp$omega, lotri({ @@ -966,6 +976,9 @@ test_that("ini(diag) and ini(-cov()) tests", { lcl ~ c(0, 0.1, 0.01, 1) })) + expect_error(mod2 %>% ini(diag(matt)), + "matt") + # Will reorder tmp <- mod2 %>% ini(diag(lcl, lvc)) expect_equal(tmp$omega, From 44c51309a1fa3510d66d93701fa893a1dc61d986 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 22:15:03 -0500 Subject: [PATCH 10/11] Fix introduced errors --- R/piping.R | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/R/piping.R b/R/piping.R index 8c3beaddd..7b9319e61 100644 --- a/R/piping.R +++ b/R/piping.R @@ -373,15 +373,16 @@ # Capture empty arguments (rxode2#688) warning("empty argument ignored") return(NULL) - } else if (is.symbol(.quoted) && - identical(.quoted, quote(`diag`))) { - .quoted <- str2lang("~diag()") + } else if (length(.quoted) == 1) { + if (identical(.quoted, quote(`diag`))) { + .quoted <- str2lang("~diag()") + } else { + .bracket[i] <- TRUE + assign(".bracket", .bracket, envir=.env) + } } else if (length(.quoted) >= 1 && identical(.quoted[[1]], quote(`diag`))) { .quoted <- as.call(c(list(quote(`~`)), .quoted)) - } else if (length(.quoted) == 1) { - .bracket[i] <- TRUE - assign(".bracket", .bracket, envir=.env) } else if (identical(.quoted[[1]], quote(`diag`))) { } else if (identical(.quoted[[1]], quote(`{`)) || From ff8e9b0ab3a86ae1b4ad19f72f1d96b1d7c0d494 Mon Sep 17 00:00:00 2001 From: Matthew Fidler <matthew.fidler@gmail.com> Date: Sat, 7 Sep 2024 22:25:55 -0500 Subject: [PATCH 11/11] Expand/fix diag() and -cov ini doc --- R/ui.R | 10 ++++++---- man/ini.Rd | 26 ++++++++++++++++++-------- 2 files changed, 24 insertions(+), 12 deletions(-) diff --git a/R/ui.R b/R/ui.R index d77738e47..5c838daad 100644 --- a/R/ui.R +++ b/R/ui.R @@ -177,10 +177,12 @@ #' label("Typical Value of Clearance (L/hr)")}. #' #' Off diagonal values of 'omega' can be set to zero using the -#' \code{diag()} or \code{diag(iivKa, iivCL)} for example removing all -#' off-diagonals can be removed with `ini(diag())`, or the off -#' diagonals between clearance and ka could be removed by -#' \code{ini(diag(iivKa, iivCL)}. +#' \code{diag()} to remove all off-diagonals can be removed with +#' `ini(diag())`. To remove covariances of 'omega' item with `iivKa`, +#' you can use `%>% ini(diag(iivKa))`. Or to remove covariances that +#' contain either `iivKa` or `iivCl` you can use `%>% ini(diag(iivKa, +#' iivCl))`. For finer control you can remove the covariance between +#' two items (like `iivKa` and `iivCl`) by `%>% ini(-cov(iivKa, iivCl)) #' #' \code{rxode2}/\code{nlmixr2} will attempt to determine some #' back-transformations for the user. For example, \code{CL <- exp(tvCL)} will diff --git a/man/ini.Rd b/man/ini.Rd index caf686237..e8dde878b 100644 --- a/man/ini.Rd +++ b/man/ini.Rd @@ -63,14 +63,17 @@ warning. 'omega' values can be set as a single value or as the values of a lower-triangular matrix. The values may be set as either a -variance-covariance matrix (the default) or as a correlation matrix for the -off-diagonals with the standard deviations on the diagonals. Names may be -set on the left side of the \code{~}. To set a variance-covariance matrix -with variance values of 2 and 3 and a covariance of -2.5 use \code{~c(2, 2.5, -3)}. To set the same matrix with names of \code{iivKa} and \code{iivCL}, use -\code{iivKa + iivCL~c(2, 2.5, 3)}. To set a correlation matrix with standard -deviations on the diagonal, use \code{cor()} like \code{iivKa + iivCL~cor(2, --0.5, 3)}. +variance-covariance matrix (the default) or as a correlation matrix +for the off-diagonals with the standard deviations on the +diagonals. Names may be set on the left side of the \code{~}. To +set a variance-covariance matrix with variance values of 2 and 3 +and a covariance of -2.5 use \code{~c(2, 2.5, 3)}. To set the same +matrix with names of \code{iivKa} and \code{iivCL}, use \code{iivKa ++ iivCL~c(2, 2.5, 3)}. To set a correlation matrix with standard +deviations on the diagonal, use \code{cor()} like \code{iivKa + +iivCL~cor(2, -0.5, 3)}. As of rxode2 3.0 you can also use +\code{iivKa ~ 2, iivCL ~ c(2.5, 3)} for covariance matrices as +well. Values may be fixed (and therefore not estimated) using either the name \code{fixed} at the end of the assignment or by calling \code{fixed()} as a @@ -99,6 +102,13 @@ estimation. The typical way to set a label so that the parameter \code{tvCL} has a label of "Typical Value of Clearance (L/hr)" is \code{tvCL <- 1; label("Typical Value of Clearance (L/hr)")}. +Off diagonal values of 'omega' can be set to zero using the +\code{diag()} to remove all off-diagonals can be removed with +\code{ini(diag())}. To remove covariances of 'omega' item with \code{iivKa}, +you can use \verb{\%>\% ini(diag(iivKa))}. Or to remove covariances that +contain either \code{iivKa} or \code{iivCl} you can use \verb{\%>\% ini(diag(iivKa, iivCl))}. For finer control you can remove the covariance between +two items (like \code{iivKa} and \code{iivCl}) by `\%>\% ini(-cov(iivKa, iivCl)) + \code{rxode2}/\code{nlmixr2} will attempt to determine some back-transformations for the user. For example, \code{CL <- exp(tvCL)} will detect that \code{tvCL} must be back-transformed by \code{exp()} for easier