Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

619 allow initial parameters to be changed to covariate pop par or between subject variability with piping #621

Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
124 changes: 124 additions & 0 deletions R/piping-ini.R
Original file line number Diff line number Diff line change
Expand Up @@ -409,6 +409,126 @@
invisible()
}

.iniHandleRecalc <- function(rxui) {
.fun <- rxUiDecompress(rxui$fun())
for (.i in ls(.fun, all=TRUE)) {
if (.i != "meta") {
assign(.i, get(.i, envir=.fun), envir=rxui)
}
}
invisible()
}

#' Handle switching theta to eta and vice versa
#'
#' This is coded as model |> ini(~par)
#'
#' @param expr Expression, this would be the ~par expression
#' @param rxui rxui uncompressed environment
#' @param envir Environment for evaluation (if needed)
#' @return Nothing, called for side effects
#' @noRd
#' @author Matthew L. Fidler
.iniHandleSwitchType <- function(expr, rxui, envir=parent.frame()) {
.var <- as.character(expr[[2]])
.iniDf <- rxui$iniDf
.w <- which(.iniDf$name == .var)
if (length(.w) != 1L) stop("cannot switch parameter type for '", .var, "'", call.=FALSE)
.theta <- .iniDf[!is.na(.iniDf$ntheta),, drop = FALSE]
.eta <- .iniDf[is.na(.iniDf$ntheta),, drop = FALSE]
if (is.na(.iniDf$ntheta[.w])) {
# switch eta to theta
.neta <- .iniDf$neta1[.w]
.eta <- .eta[.eta$neta1 != .neta,, drop = FALSE]
.eta <- .eta[.eta$neta2 != .neta,, drop = FALSE]
.eta$neta1 <- .eta$neta1 - ifelse(.eta$neta1 < .neta, 0L, 1L)
.eta$neta2 <- .eta$neta2 - ifelse(.eta$neta2 < .neta, 0L, 1L)
.newTheta <- .iniDf[.w, ]
.newTheta$neta1 <- NA_integer_
.newTheta$neta2 <- NA_integer_
if (length(.theta$ntheta) == 0L) {
.newTheta$ntheta <- 1L
} else {
.newTheta$ntheta <- max(.theta$ntheta) + 1L
}
.minfo(paste0("convert '", .var, "' from between subject variability to population parameter"))
.theta <- rbind(.theta, .newTheta)
} else {
# switch theta to eta
if (!is.na(.iniDf$err[.w])) {
stop("cannot switch error parameter '", .var,
"' to a different type", call. = FALSE)
}
.ntheta <- .iniDf$ntheta[.w]
.theta <- .theta[.theta$ntheta != .ntheta,, drop = FALSE]
.theta$ntheta <- .theta$ntheta - ifelse(.theta$ntheta < .ntheta, 0L, 1L)
.newEta <- .iniDf[.w, ]
.newEta$ntheta <- NA_integer_
if (length(.eta$neta1) == 0L) {
.newEta$neta1 <- .newEta$neta2 <- 1L
} else {
.newEta$neta1 <- .newEta$neta2 <- max(.eta$neta1) + 1L

Check warning on line 470 in R/piping-ini.R

View check run for this annotation

Codecov / codecov/patch

R/piping-ini.R#L470

Added line #L470 was not covered by tests
}
.minfo(paste0("convert '", .var, "' from population parameter to between subject variability"))
if (.newEta$est == 0) {
.minfo("old initial estimate is zero, changing to 1")
.newEta$est <- 1
} else if (.newEta$est < 0) {
.minfo("old initial estimate was negative, changing to positive")
.newEta$est <- -.newEta$est
}
.newEta$lower <- -Inf
.newEta$upper <- Inf
.eta <- rbind(.eta, .newEta)
}
.ini <- rbind(.theta, .eta)
assign("iniDf", .ini, envir=rxui)
# This may change mu-referencing, recalculate everything
.iniHandleRecalc(rxui)
invisible()
}

#' Handle dropping parameter and treating as if it is a covariate
#'
#' This is coded as model |> ini(-par)
#'
#' @param expr Expression, this would be the ~par expression
#' @param rxui rxui uncompressed environment
#' @param envir Environment for evaluation (if needed)
#' @return Nothing, called for side effects
#' @noRd
#' @author Matthew L. Fidler
.iniHandleDropType <- function(expr, rxui, envir=parent.frame()) {
.var <- as.character(expr[[2]])
.iniDf <- rxui$iniDf
.w <- which(.iniDf$name == .var)
if (length(.w) != 1L) stop("no initial estimates for '", .var, "', cannot change to covariate", call.=FALSE)

Check warning on line 505 in R/piping-ini.R

View check run for this annotation

Codecov / codecov/patch

R/piping-ini.R#L505

Added line #L505 was not covered by tests
.theta <- .iniDf[!is.na(.iniDf$ntheta),, drop = FALSE]
.eta <- .iniDf[is.na(.iniDf$ntheta),, drop = FALSE]
if (is.na(.iniDf$ntheta[.w])) {
.minfo(paste0("changing between subject variability parameter '", .var, "' to covariate parameter"))
.neta <- .iniDf$neta1[.w]
.eta <- .eta[.eta$neta1 != .neta,, drop = FALSE]
.eta <- .eta[.eta$neta2 != .neta,, drop = FALSE]
.eta$neta1 <- .eta$neta1 - ifelse(.eta$neta1 < .neta, 0L, 1L)
.eta$neta2 <- .eta$neta2 - ifelse(.eta$neta2 < .neta, 0L, 1L)
} else {
if (!is.na(.iniDf$err[.w])) {
stop("cannot switch error parameter '", .var,
"' to a covariate", call. = FALSE)

Check warning on line 518 in R/piping-ini.R

View check run for this annotation

Codecov / codecov/patch

R/piping-ini.R#L516-L518

Added lines #L516 - L518 were not covered by tests
}
.minfo(paste0("changing population parameter '", .var, "' to covariate parameter"))
.ntheta <- .iniDf$ntheta[.w]
.theta <- .theta[.theta$ntheta != .ntheta,, drop = FALSE]
.theta$ntheta <- .theta$ntheta - ifelse(.theta$ntheta < .ntheta, 0L, 1L)

Check warning on line 523 in R/piping-ini.R

View check run for this annotation

Codecov / codecov/patch

R/piping-ini.R#L520-L523

Added lines #L520 - L523 were not covered by tests
}
.ini <- rbind(.theta, .eta)
assign("iniDf", .ini, envir=rxui)
# This will change covariates, recalculate everything
.iniHandleRecalc(rxui)
invisible()
}

#' Update the iniDf of a model
#'
#' @param expr Expression for parsing
Expand Down Expand Up @@ -464,6 +584,10 @@
expr[[3]] <- eval(as.call(list(quote(`lotri`), as.call(list(quote(`{`), expr)))),
envir=envir)[1, 1]
.iniHandleFixOrUnfixEqual(expr=expr, rxui=rxui, envir=envir, maxLen=1L)
} else if (.isTildeExpr(expr)) {
.iniHandleSwitchType(expr=expr, rxui=rxui, envir=envir)
} else if (.isIniDropExpression(expr)) {
.iniHandleDropType(expr=expr, rxui=rxui, envir=envir)
} else {
# Can this error be improved to clarify what is the expression causing the
# issue? It needs a single character string representation of something
Expand Down
8 changes: 8 additions & 0 deletions R/piping-model.R
Original file line number Diff line number Diff line change
Expand Up @@ -180,6 +180,14 @@ model.rxModelVars <- model.rxode2
.matchesLangTemplate(expr, str2lang(". = ."))
}

.isTildeExpr <- function(expr) {
.matchesLangTemplate(expr, str2lang("~ ."))
}

.isIniDropExpression <- function(expr) {
.matchesLangTemplate(expr, str2lang("- ."))
}

# get the left hand side of an assignment or endpoint; returns NULL if the input
# is not an assignment or endpoint
.getLhs <- function(expr) {
Expand Down
10 changes: 6 additions & 4 deletions R/piping.R
Original file line number Diff line number Diff line change
Expand Up @@ -295,10 +295,12 @@
} else if (identical(.quoted[[1]], quote(`as.formula`))) {
.quoted <- .quoted[[2]]
} else if (identical(.quoted[[1]], quote(`~`))) {
.quoted[[3]] <- .iniSimplifyFixUnfix(.quoted[[3]])
if (identical(.quoted[[3]], quote(`fix`)) ||
identical(.quoted[[3]], quote(`unfix`))) {
.quoted <- as.call(list(quote(`<-`), .quoted[[2]], .quoted[[3]]))
if (length(.quoted) == 3L) {
.quoted[[3]] <- .iniSimplifyFixUnfix(.quoted[[3]])
if (identical(.quoted[[3]], quote(`fix`)) ||
identical(.quoted[[3]], quote(`unfix`))) {
.quoted <- as.call(list(quote(`<-`), .quoted[[2]], .quoted[[3]]))
}
}
} else if (identical(.quoted[[1]], quote(`$`))) {
.tmp <- try(eval(.quoted), silent=TRUE)
Expand Down
168 changes: 168 additions & 0 deletions tests/testthat/test-piping-ini.R
Original file line number Diff line number Diff line change
Expand Up @@ -569,3 +569,171 @@ test_that("append allows promoting from covariate (#472)", {
)
expect_equal(newmod$iniDf$name, c("lka", "ka_dose", "lcl", "lvc", "propSd"))
})

test_that("change ini type with ~", {

mod <- function() {
ini({
lka <- 0.45
lcl <- 1
lvc <- 3.45
propSd <- 0.5
})
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
cp ~ prop(propSd)
})
}

mod1 <- mod |> ini( ~ lka)
expect_equal(mod1$omega, lotri(lka ~ 0.45))

mod2 <- mod1 |> ini( ~ lka)
expect_equal(mod2$omega, NULL)

expect_error(mod1 |> ini( ~ propSd))

expect_error(mod1 |> ini( ~ matt))

## all etas

mod <- function() {
ini({
lka ~ 0.45
lcl ~ 1
lvc ~ 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
})
}

mod2 <- mod |> ini( ~ lka)

expect_equal(mod2$omega, lotri(lcl ~ 1, lvc ~ 3.45))

# remove correlated eta

mod <- function() {
ini({
lka + lcl + lvc ~
c(0.45,
0.01, 1,
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
})
}

mod2 <- mod |> ini( ~ lka)

expect_equal(mod2$omega, lotri(lcl + lvc ~ c(1,
-0.01, 3.45)))


# negative and zero

mod <- function() {
ini({
lka <- 0.45
lcl <- -1
lvc <- 0
})
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
})
}

mod2 <- mod |> ini( ~ lcl)

expect_equal(mod2$omega, lotri(lcl ~ 1))

mod2 <- mod |> ini( ~ lvc)

expect_equal(mod2$omega, lotri(lvc ~ 1))

mod3 <- mod2 |> ini( ~ lvc)

expect_equal(mod3$omega, NULL)

mod4 <- mod3 |> ini( ~ lvc)

expect_equal(mod4$omega, lotri(lvc ~ 1))

})



test_that("change ini variable to covariate with -", {

mod <- function() {
ini({
lka + lcl + lvc ~
c(0.45,
0.01, 1,
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
})
}

mod2 <- mod |> ini(-lka)

expect_equal(mod2$allCovs, "lka")
expect_equal(mod2$omega, lotri(lcl + lvc ~ c(1, -0.01, 3.45)))

mod <- function() {
ini({
lka ~ 0.45
lcl ~ 1
lvc ~ 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
})
}

mod2 <- mod |> ini(-lka)

expect_equal(mod2$allCovs, "lka")


})
Loading