Closed DrOJG closed 2 months ago
billdenney
commented
4 months ago I don't think that there are any known version incompatibilities. But, changing libraries may possibly require recompiling the nlmixr2 tools because of linking changes. Can you try reinstalling nlmixr2 and see if it works then?
DrOJG
commented
4 months ago Hi, After uninstalling nlmixr2 and all dependencies, purging the renv cache, and reinstalling, the problem is still persisting. Happy to provide any info that might help with solving. Best, O.
mattfidler
commented
4 months ago This is likely not the case for all models. If there is an error in the table step, nlmixr will return an object without the table step.
If the data are not propritary and the model is not propriatary the easiest way to provide a reproducible example is to save the model and then use that to try to add tables again. I believe the function is calcTable() but I dont currently have access to a computer so this is based on memory.
Once I have a reproducible example I can try to figure out what is going on in your case.
DrOJG
commented
4 months ago Hi Matt, At this point I would prefer not to share the dataset or the complete model as we are looking to publish soon, I apologise for this. I have however done some further investigation which will hopefully be useful in diagnosing.
I am currently running the model fitting with the below command:
model.fit <- nlmixr2(model, data_set, est = "foce", table = tableControl(cwres = TRUE, npde = TRUE))As stated previously, this is returning a fit object without the table step. When I try to add a table to the object with the addTable(), it throws the following error:
→ Calculating residuals/tables
[====|====|====|====|====|====|====|====|====|====] 0:00:00
Error in `$<-.data.frame`(`*tmp*`, "sim.id", value = numeric(0)) :
replacement has 0 rows, data has 302
In addition: Warning message:
In max(.sim$sim.id) : no non-missing arguments to max; returning -InfAdditionally, I have tried running both the Friberg myelosuppression example model and the Nimotuzumab example model in the same environment and both returned model fit dataframes. It is clearly a specific issue, but strange that this code previously ran fine with no alterations to the model or dataset.
Thanks, O.
mattfidler
commented
4 months ago If you can add a traceback it would be more helpful. Currently I have insufficient information to do anything about this ☹️
mattfidler
commented
4 months ago This is done with traceback() from the place after addTable(). This way it gives where the error occurs without providing the dataset or model.
DrOJG
commented
4 months ago Sorry Matt, I have been on leave so was not able to respond! The following command:
model_with_table <- addTable(model.fit.object, table = tableControl(cwres = TRUE, npde = TRUE))throws the aforementioned error, the traceback() is as follows:
> traceback()
11: stop(sprintf(ngettext(N, "replacement has %d row, data has %d",
"replacement has %d rows, data has %d"), N, nrows), domain = NA)
10: `$<-.data.frame`(`*tmp*`, "sim.id", value = numeric(0))
9: `$<-`(`*tmp*`, "sim.id", value = numeric(0))
8: vpcSim(fit, n = table$nsim, seed = table$seed, addDosing = addDosing,
subsetNonmem = subsetNonmem)
7: .calcCwres0(fit, data, thetaEtaParameters, table, dv = dv, predOnly,
addDosing, subsetNonmem, keep, npde, .prdLst = .prdLst)
6: .calcCwres(fit, data = fit$dataSav, thetaEtaParameters = .thetaEtaParameters,
table = table, dv = .ret[[1]][[1]], predOnly = .predOnly,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
keep = keep, .prdLst = .prdLst, npde = .npde2)
5: .calcTables(.fit, data = data, table = table, keep = keep)
4: force(code)
3: nlmixrWithTiming("table", {
keep <- unique(c(keep, "nlmixrRowNums"))
.malert("Calculating residuals/tables")
.objName <- substitute(object)
if (!inherits(object, "nlmixr2FitCore")) {
stop("requires a nlmixr2 fit object", call. = FALSE)
}
.fit <- object$env
if (exists("origControl", .fit)) {
.control <- .fit$origControl
}
else if (exists("control", .fit)) {
.control <- .fit$control
}
else {
.control <- foceiControl()
}
if (is.null(.fit$omega)) {
.df <- .calcIres(.fit, data = data, table = table, dv = NULL,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
...
2: nlmixr2est::addTable(object = object, updateObject = updateObject,
data = data, thetaEtaParameters = thetaEtaParameters, table = table,
keep = keep, drop = drop, envir = envir)
1: addTable(model.fit.object, table = tableControl(cwres = TRUE,
npde = TRUE))Let me know if this gives any insight. I am also conducting my own investigations into this to narrow down what specifically is causing the error.
p.s. Thank you for your efforts here, I appreaciate it!
DrOJG
commented
4 months ago Hi Matt, To update further, I appear to have isolated the issue.
A single covariate, coded in the form tL0 + beta_pbcb_L0 * PB_or_CB + eta.L0 where PB_or_CB is a binary indicator variable in the dataset, was causing the issue. Removing specifically this covariate from the model fixed the issue, albeit changing the fit slightly.
I am still not sure why this was causing an issue, as it was previously included without causing problems, but it does appear to be the source of the error. It may be it was causing some instability in the fit but I am not sure what changed to cause it to break the table calculation.
Thanks, O.
billdenney
commented
4 months ago That is very helpful diagnosis. Is the covariate column always or almost always the same value?
DrOJG
commented
4 months ago Hi Bill,
As I say, it is a simple binary covariate (0 or 1), but I have just checked through the dataset and found it is set to 1 in approximately 20% of rows.
The parameter estimate when included was approx. -0.259, but has a fairly wide confidence interval of (-0.481, -0.0373)
Thanks,
O.
mattfidler
commented
4 months ago Hi @DrOJG
I believe that this happens with the vpcSim() function that simulates both the vpc and the information for npde. As a confirmation you remove npde it is likely that this will return a data frame.
The problem is with this covariate (and its large range), sometimes it selects models that cannot be solved. In general, these simulations are rejected and then another simulation replaces the simulations that could not be solved.
You have reached an edge case where this replacement isn't working well. I will try to help you debug this further. First I will transfer this to nlmixr2est where the vpcSim() model lives.
mattfidler
commented
4 months ago Just before you model_with_table <- addTable(model.fit.object, table = tableControl(cwres = TRUE, npde = TRUE)), add the following code:
vpcSim <- function(object, ..., keep=NULL, n=300,
pred=FALSE, seed=1009, nretry=50, minN=10,
normRelated=TRUE) {
checkmate::assertIntegerish(minN, len=1, any.missing=FALSE, lower=2)
checkmate::assertLogical(pred, len=1, any.missing=FALSE)
checkmate::assertIntegerish(nretry, len=1, any.missing=FALSE, lower=0)
checkmate::assertLogical(normRelated, len=1, any.missing=FALSE)
checkmate::assertCharacter(keep, null.ok=TRUE, pattern="^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$")
checkmate::assertIntegerish(seed)
assignInMyNamespace(".finalUiCompressed", FALSE)
on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
set.seed(seed)
.si <- object$simInfo
.si$object <- eval(.getSimModel(object, hideIpred=FALSE))
.w <- which(names(.si) == "rx")
.si <- .si[-.w]
.si$nsim <- n
.si <- c(.si, list(...))
.pt <- proc.time()
.si$keep <- unique(c(keep, "nlmixrRowNums"))
.data <- .si$events
.data$nlmixrRowNums <- seq_along(.data[, 1])
if (normRelated) {
.ui <- object$ui
.predDf <-.ui$predDf
if (all(.predDf$dist %in% c("norm", "dnorm","t", "cauchy"))) {
} else {
if (is.null(.data$CMT)) {
.ds <- object$dataSav
.data$nlmixrRowNums <- seq_along(.data[,1])
.ds <- .ds[, c("CMT", "nlmixrRowNums")]
.data <- merge(.data, .ds, by ="nlmixrRowNums")
.data <- .data[order(.data$nlmixrRowNums),]
}
.lst <- .Call(`_nlmixr2est_filterNormalLikeAndDoses`,
.data$CMT, .predDf$distribution, .predDf$cmt)
.lst$nlmixrRowNums <- .data[.lst$filter, "nlmixrRowNums"]
if (.lst$nnorm == 0L) {
stop("need normal data for vpcSim (or use normRelated=FALSE)")
}
.data <- .data[.lst$filter, ]
}
}
.si$events <- .data
.si$thetaMat <- NULL
.si$dfSub <- NULL
.si$dfObs <- NULL
.si$returnType <- "data.frame.TBS"
.sim <- do.call(rxode2::rxSolve, .si)
# now look for how many have missing values
.w <- which(is.na(.sim$ipred))
.nretry <- 0
while (length(.w) > 0 && .nretry < nretry) {
.w <- which(is.na(.sim$ipred))
.simIds <- unique(.sim$sim.id[.w])
.sim <- .sim[!(.sim$sim.id %in% .simIds), ]
print("loc1:")
print(as.integer(factor(.sim$sim.id)))
.sim$sim.id <- as.integer(factor(.sim$sim.id))
.mx <- max(.sim$sim.id)
.si$nsim <- n - .mx
.adjust <- FALSE
if (.si$nsim < minN) {
.si$nsim <- minN
.adjust <- TRUE
}
.sim2 <- do.call(rxode2::rxSolve, .si)
if (!("sim.id" %in% names(.sim2))) {
.sim2$sim.id <- .simIds[1]
}
if (.adjust) {
# Select simulations without NA ipreds in them
.w <- which(is.na(.sim2$ipred))
.simIds <- unique(.sim2$sim.id[.w])
.allSimIds <- unique(.sim2$sim.id)
.simIds <- .allSimIds[!(.allSimIds %in% .simIds)]
.simIds <- .simIds[seq_len(min(n - .mx, length(.simIds)))]
.sim2 <- .sim2[.sim2$sim.id %in% .simIds, ]
}
print("loc2:")
print(.sim2$sim.id + .mx)
.sim2$sim.id <- .sim2$sim.id + .mx
.sim <- rbind(.sim, .sim2)
.w <- which(is.na(.sim$ipred))
.nretry <- .nretry + 1
}
if (.nretry != 0) {
if (length(.w) == 0) {
warning("'NA' values in vpc or npde simulation, re-simulated until all simulations were successful",
call.=FALSE)
} else {
warning("'NA' values in vpc or npde simulation",
call.=FALSE)
}
}
if (pred) {
.si$nsim <- n # restore for pred
.si2 <- .si
.si2$params <- c(
.si$params, setNames(rep(0, dim(.si$omega)[1]), dimnames(.si$omega)[[2]]),
setNames(rep(0, dim(.si$sigma)[1]), dimnames(.si$sigma)[[2]]))
.si2$omega <- NULL
.si2$sigma <- NULL
.si2$returnType <- "data.frame"
.si2$nStud <- 1
.si2$nsim <- NULL
assignInMyNamespace(".lastPredSimulationInfo", .si2)
.sim2 <- do.call(rxode2::rxSolve, .si2)
.sim$pred <- .sim2$sim
}
.sim <- vpcNameDataCmts(object, .sim)
.cls <- c("nlmixr2vpcSim", class(.sim))
.fit <- object
.cls0 <- c("rxHidden", class(.fit))
attr(.cls0, ".foceiEnv") <- attr(class(.fit), ".foceiEnv")
class(.fit) <- .cls0
attr(.cls, "fit") <- .fit
class(.sim) <- .cls
return(.sim)
}
environment(vpcSim) <- asNamespace('nlmixr2est')
assignInNamespace("vpcSim", vpcSim, ns = "nlmixr2est")This may fix your problem, though if not it will print out some information that will help investigate the issue further.
mattfidler
commented
4 months ago Let me know if it fixes it, as if it does, I will add this to nlmixr2est to fix this in future revisions.
DrOJG
commented
4 months ago Apologies in advance for the long comment!
As a confirmation you remove npde it is likely that this will return a data frame.
Yes, this does indeed return a model fit object with the data frame.
Just before you model_with_table <- addTable(model.fit.object, table = tableControl(cwres = TRUE, npde = TRUE)), add the following code...
So addTable() after running the code you provided game me the following:
And so on, same warning x50.
traceback() gives me:
> traceback()
7: .calcCwres0(fit, data, thetaEtaParameters, table, dv = dv, predOnly,
addDosing, subsetNonmem, keep, npde, .prdLst = .prdLst)
6: .calcCwres(fit, data = fit$dataSav, thetaEtaParameters = .thetaEtaParameters,
table = table, dv = .ret[[1]][[1]], predOnly = .predOnly,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
keep = keep, .prdLst = .prdLst, npde = .npde2)
5: .calcTables(.fit, data = data, table = table, keep = keep)
4: force(code)
3: nlmixrWithTiming("table", {
keep <- unique(c(keep, "nlmixrRowNums"))
.malert("Calculating residuals/tables")
.objName <- substitute(object)
if (!inherits(object, "nlmixr2FitCore")) {
stop("requires a nlmixr2 fit object", call. = FALSE)
}
.fit <- object$env
if (exists("origControl", .fit)) {
.control <- .fit$origControl
}
else if (exists("control", .fit)) {
.control <- .fit$control
}
else {
.control <- foceiControl()
}
if (is.null(.fit$omega)) {
.df <- .calcIres(.fit, data = data, table = table, dv = NULL,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
...
2: nlmixr2est::addTable(object = object, updateObject = updateObject,
data = data, thetaEtaParameters = thetaEtaParameters, table = table,
keep = keep, drop = drop, envir = envir)
1: addTable(vpc_tshoot_1, table = tableControl(cwres = TRUE, npde = TRUE))Hope this is helpful!
mattfidler
commented
4 months ago Ok, try the following instead:
vpcSim <- function(object, ..., keep=NULL, n=300,
pred=FALSE, seed=1009, nretry=50, minN=10,
normRelated=TRUE) {
checkmate::assertIntegerish(minN, len=1, any.missing=FALSE, lower=2)
checkmate::assertLogical(pred, len=1, any.missing=FALSE)
checkmate::assertIntegerish(nretry, len=1, any.missing=FALSE, lower=0)
checkmate::assertLogical(normRelated, len=1, any.missing=FALSE)
checkmate::assertCharacter(keep, null.ok=TRUE, pattern="^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$")
checkmate::assertIntegerish(seed)
assignInMyNamespace(".finalUiCompressed", FALSE)
on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
set.seed(seed)
.si <- object$simInfo
.si$object <- eval(.getSimModel(object, hideIpred=FALSE))
.w <- which(names(.si) == "rx")
.si <- .si[-.w]
.si$nsim <- n
.si <- c(.si, list(...))
.pt <- proc.time()
.si$keep <- unique(c(keep, "nlmixrRowNums"))
.data <- .si$events
.data$nlmixrRowNums <- seq_along(.data[, 1])
if (normRelated) {
.ui <- object$ui
.predDf <-.ui$predDf
if (all(.predDf$dist %in% c("norm", "dnorm","t", "cauchy"))) {
} else {
if (is.null(.data$CMT)) {
.ds <- object$dataSav
.data$nlmixrRowNums <- seq_along(.data[,1])
.ds <- .ds[, c("CMT", "nlmixrRowNums")]
.data <- merge(.data, .ds, by ="nlmixrRowNums")
.data <- .data[order(.data$nlmixrRowNums),]
}
.lst <- .Call(`_nlmixr2est_filterNormalLikeAndDoses`,
.data$CMT, .predDf$distribution, .predDf$cmt)
.lst$nlmixrRowNums <- .data[.lst$filter, "nlmixrRowNums"]
if (.lst$nnorm == 0L) {
stop("need normal data for vpcSim (or use normRelated=FALSE)")
}
.data <- .data[.lst$filter, ]
}
}
.si$events <- .data
.si$thetaMat <- NULL
.si$dfSub <- NULL
.si$dfObs <- NULL
.si$returnType <- "data.frame.TBS"
.sim <- do.call(rxode2::rxSolve, .si)
# now look for how many have missing values
.w <- which(is.na(.sim$ipred))
.nretry <- 0
while (length(.w) > 0 && .nretry < nretry) {
.w <- which(is.na(.sim$ipred))
.simIds <- unique(.sim$sim.id[.w])
.sim <- .sim[!(.sim$sim.id %in% .simIds),, drop = FALSE]
if (length(.sim$sim.id) == 0) {
if (length(.w) > 0) {
warning("some NA ipreds replaced with zero; couldn't figure out how to re-simulate",
call.=FALSE)
.sim$ipred[.w] <- 0
}
break
}
.sim$sim.id <- as.integer(factor(.sim$sim.id))
.mx <- max(.sim$sim.id)
.si$nsim <- n - .mx
.adjust <- FALSE
if (.si$nsim < minN) {
.si$nsim <- minN
.adjust <- TRUE
}
.sim2 <- do.call(rxode2::rxSolve, .si)
if (!("sim.id" %in% names(.sim2))) {
.sim2$sim.id <- .simIds[1]
}
if (.adjust) {
# Select simulations without NA ipreds in them
.w <- which(is.na(.sim2$ipred))
.simIds <- unique(.sim2$sim.id[.w])
.allSimIds <- unique(.sim2$sim.id)
.simIds <- .allSimIds[!(.allSimIds %in% .simIds)]
.simIds <- .simIds[seq_len(min(n - .mx, length(.simIds)))]
.sim2 <- .sim2[.sim2$sim.id %in% .simIds, ]
}
.sim2$sim.id <- .sim2$sim.id + .mx
.sim <- rbind(.sim, .sim2)
.w <- which(is.na(.sim$ipred))
.nretry <- .nretry + 1
}
if (.nretry != 0) {
if (length(.w) == 0) {
warning("'NA' values in vpc or npde simulation, re-simulated until all simulations were successful",
call.=FALSE)
} else {
warning("'NA' values in vpc or npde simulation",
call.=FALSE)
}
}
if (pred) {
.si$nsim <- n # restore for pred
.si2 <- .si
.si2$params <- c(
.si$params, setNames(rep(0, dim(.si$omega)[1]), dimnames(.si$omega)[[2]]),
setNames(rep(0, dim(.si$sigma)[1]), dimnames(.si$sigma)[[2]]))
.si2$omega <- NULL
.si2$sigma <- NULL
.si2$returnType <- "data.frame"
.si2$nStud <- 1
.si2$nsim <- NULL
assignInMyNamespace(".lastPredSimulationInfo", .si2)
.sim2 <- do.call(rxode2::rxSolve, .si2)
.sim$pred <- .sim2$sim
}
.sim <- vpcNameDataCmts(object, .sim)
.cls <- c("nlmixr2vpcSim", class(.sim))
.fit <- object
.cls0 <- c("rxHidden", class(.fit))
attr(.cls0, ".foceiEnv") <- attr(class(.fit), ".foceiEnv")
class(.fit) <- .cls0
attr(.cls, "fit") <- .fit
class(.sim) <- .cls
return(.sim)
}
environment(vpcSim) <- asNamespace('nlmixr2est')
assignInNamespace("vpcSim", vpcSim, ns = "nlmixr2est")
That is use this vpcSim() before you calculate your tables. In theory, this should fix the problem, but let me know.
DrOJG
commented
4 months ago So after running the new vpcSim() code above, addTable() gives the following error:
> mod_with_table <- addTable(vpc_tshoot_1, table = tableControl(cwres = TRUE, npde = TRUE))
→ Calculating residuals/tables
[====|====|====|====|====|====|====|====|====|====] 0:00:00
Error in `$<-.data.frame`(`*tmp*`, ipred, value = c(0, 0, 0, 0, 0, 0, :
replacement has 90600 rows, data has 0
In addition: Warning message:
some NA ipreds replaced with zero; couldn't figure out how to re-simulate traceback():
> traceback()
11: stop(sprintf(ngettext(N, "replacement has %d row, data has %d",
"replacement has %d rows, data has %d"), N, nrows), domain = NA)
10: `$<-.data.frame`(`*tmp*`, ipred, value = c(0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
... at #61
9: `$<-`(`*tmp*`, ipred, value = c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
... at #61
8: vpcSim(fit, n = table$nsim, seed = table$seed, addDosing = addDosing,
subsetNonmem = subsetNonmem)
7: .calcCwres0(fit, data, thetaEtaParameters, table, dv = dv, predOnly,
addDosing, subsetNonmem, keep, npde, .prdLst = .prdLst)
6: .calcCwres(fit, data = fit$dataSav, thetaEtaParameters = .thetaEtaParameters,
table = table, dv = .ret[[1]][[1]], predOnly = .predOnly,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
keep = keep, .prdLst = .prdLst, npde = .npde2)
5: .calcTables(.fit, data = data, table = table, keep = keep)
4: force(code)
3: nlmixrWithTiming("table", {
keep <- unique(c(keep, "nlmixrRowNums"))
.malert("Calculating residuals/tables")
.objName <- substitute(object)
if (!inherits(object, "nlmixr2FitCore")) {
stop("requires a nlmixr2 fit object", call. = FALSE)
}
.fit <- object$env
if (exists("origControl", .fit)) {
.control <- .fit$origControl
}
else if (exists("control", .fit)) {
.control <- .fit$control
}
else {
.control <- foceiControl()
}
if (is.null(.fit$omega)) {
.df <- .calcIres(.fit, data = data, table = table, dv = NULL,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
...
2: nlmixr2est::addTable(object = object, updateObject = updateObject,
data = data, thetaEtaParameters = thetaEtaParameters, table = table,
keep = keep, drop = drop, envir = envir)
1: addTable(vpc_tshoot_1, table = tableControl(cwres = TRUE, npde = TRUE))
I really don't know why it is flagging bad simulations and then is not able to find them in your problem. I'm unsure how much I can do without data.
However, perhaps this will work?
vpcSim <- function(object, ..., keep=NULL, n=300,
pred=FALSE, seed=1009, nretry=50, minN=10,
normRelated=TRUE) {
checkmate::assertIntegerish(minN, len=1, any.missing=FALSE, lower=2)
checkmate::assertLogical(pred, len=1, any.missing=FALSE)
checkmate::assertIntegerish(nretry, len=1, any.missing=FALSE, lower=0)
checkmate::assertLogical(normRelated, len=1, any.missing=FALSE)
checkmate::assertCharacter(keep, null.ok=TRUE, pattern="^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$")
checkmate::assertIntegerish(seed)
assignInMyNamespace(".finalUiCompressed", FALSE)
on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
set.seed(seed)
.si <- object$simInfo
.si$object <- eval(.getSimModel(object, hideIpred=FALSE))
.w <- which(names(.si) == "rx")
.si <- .si[-.w]
.si$nsim <- n
.si <- c(.si, list(...))
.pt <- proc.time()
.si$keep <- unique(c(keep, "nlmixrRowNums"))
.data <- .si$events
.data$nlmixrRowNums <- seq_along(.data[, 1])
if (normRelated) {
.ui <- object$ui
.predDf <-.ui$predDf
if (all(.predDf$dist %in% c("norm", "dnorm","t", "cauchy"))) {
} else {
if (is.null(.data$CMT)) {
.ds <- object$dataSav
.data$nlmixrRowNums <- seq_along(.data[,1])
.ds <- .ds[, c("CMT", "nlmixrRowNums")]
.data <- merge(.data, .ds, by ="nlmixrRowNums")
.data <- .data[order(.data$nlmixrRowNums),]
}
.lst <- .Call(`_nlmixr2est_filterNormalLikeAndDoses`,
.data$CMT, .predDf$distribution, .predDf$cmt)
.lst$nlmixrRowNums <- .data[.lst$filter, "nlmixrRowNums"]
if (.lst$nnorm == 0L) {
stop("need normal data for vpcSim (or use normRelated=FALSE)")
}
.data <- .data[.lst$filter, ]
}
}
.si$events <- .data
.si$thetaMat <- NULL
.si$dfSub <- NULL
.si$dfObs <- NULL
.si$returnType <- "data.frame.TBS"
.sim <- do.call(rxode2::rxSolve, .si)
if (!("sim.id" %in% names(.sim))) {
.sim2$sim.id <- 1
}
# now look for how many have missing values
.w <- which(is.na(.sim$ipred))
.nretry <- 0
while (length(.w) > 0 && .nretry < nretry) {
.w <- which(is.na(.sim$ipred))
.simIds <- unique(.sim$sim.id[.w])
.sim <- .sim[!(.sim$sim.id %in% .simIds),, drop = FALSE]
if (length(.sim$sim.id) == 0) {
warning("when filtering for simulations, could not find any though some were flagged, be cautious with results",
call.=FALSE)
break
}
.sim$sim.id <- as.integer(factor(.sim$sim.id))
.mx <- max(.sim$sim.id)
.si$nsim <- n - .mx
.adjust <- FALSE
if (.si$nsim < minN) {
.si$nsim <- minN
.adjust <- TRUE
}
.sim2 <- do.call(rxode2::rxSolve, .si)
if (!("sim.id" %in% names(.sim2))) {
.sim2$sim.id <- .simIds[1]
}
if (.adjust) {
# Select simulations without NA ipreds in them
.w <- which(is.na(.sim2$ipred))
.simIds <- unique(.sim2$sim.id[.w])
.allSimIds <- unique(.sim2$sim.id)
.simIds <- .allSimIds[!(.allSimIds %in% .simIds)]
.simIds <- .simIds[seq_len(min(n - .mx, length(.simIds)))]
.sim2 <- .sim2[.sim2$sim.id %in% .simIds, ]
}
.sim2$sim.id <- .sim2$sim.id + .mx
.sim <- rbind(.sim, .sim2)
.w <- which(is.na(.sim$ipred))
.nretry <- .nretry + 1
}
if (.nretry != 0) {
if (length(.w) == 0) {
warning("'NA' values in vpc or npde simulation, re-simulated until all simulations were successful",
call.=FALSE)
} else {
warning("'NA' values in vpc or npde simulation",
call.=FALSE)
}
}
if (pred) {
.si$nsim <- n # restore for pred
.si2 <- .si
.si2$params <- c(
.si$params, setNames(rep(0, dim(.si$omega)[1]), dimnames(.si$omega)[[2]]),
setNames(rep(0, dim(.si$sigma)[1]), dimnames(.si$sigma)[[2]]))
.si2$omega <- NULL
.si2$sigma <- NULL
.si2$returnType <- "data.frame"
.si2$nStud <- 1
.si2$nsim <- NULL
assignInMyNamespace(".lastPredSimulationInfo", .si2)
.sim2 <- do.call(rxode2::rxSolve, .si2)
.sim$pred <- .sim2$sim
}
.sim <- vpcNameDataCmts(object, .sim)
.cls <- c("nlmixr2vpcSim", class(.sim))
.fit <- object
.cls0 <- c("rxHidden", class(.fit))
attr(.cls0, ".foceiEnv") <- attr(class(.fit), ".foceiEnv")
class(.fit) <- .cls0
attr(.cls, "fit") <- .fit
class(.sim) <- .cls
return(.sim)
}
environment(vpcSim) <- asNamespace('nlmixr2est')
assignInNamespace("vpcSim", vpcSim, ns = "nlmixr2est")Now throwing this error:
> mod_with_table <- addTable(vpc_tshoot_2, table = tableControl(cwres = TRUE, npde = TRUE))
→ Calculating residuals/tables
[====|====|====|====|====|====|====|====|====|====] 0:00:00
Error: Col::subvec(): indices out of bounds or incorrectly used
In addition: Warning message:
when filtering for simulations, could not find any though some were flagged, be cautious with results
traceback() gives:
> traceback()
8: stop(structure(list(message = "Col::subvec(): indices out of bounds or incorrectly used",
call = NULL, cppstack = NULL), class = c("std::out_of_range",
"C++Error", "error", "condition")))
7: .calcCwres0(fit, data, thetaEtaParameters, table, dv = dv, predOnly,
addDosing, subsetNonmem, keep, npde, .prdLst = .prdLst)
6: .calcCwres(fit, data = fit$dataSav, thetaEtaParameters = .thetaEtaParameters,
table = table, dv = .ret[[1]][[1]], predOnly = .predOnly,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
keep = keep, .prdLst = .prdLst, npde = .npde2)
5: .calcTables(.fit, data = data, table = table, keep = keep)
4: force(code)
3: nlmixrWithTiming("table", {
keep <- unique(c(keep, "nlmixrRowNums"))
.malert("Calculating residuals/tables")
.objName <- substitute(object)
if (!inherits(object, "nlmixr2FitCore")) {
stop("requires a nlmixr2 fit object", call. = FALSE)
}
.fit <- object$env
if (exists("origControl", .fit)) {
.control <- .fit$origControl
}
else if (exists("control", .fit)) {
.control <- .fit$control
}
else {
.control <- foceiControl()
}
if (is.null(.fit$omega)) {
.df <- .calcIres(.fit, data = data, table = table, dv = NULL,
addDosing = table$addDosing, subsetNonmem = table$subsetNonmem,
...
2: nlmixr2est::addTable(object = object, updateObject = updateObject,
data = data, thetaEtaParameters = thetaEtaParameters, table = table,
keep = keep, drop = drop, envir = envir)
1: addTable(vpc_tshoot_2, table = tableControl(cwres = TRUE, npde = TRUE))Out of interest, is there an easy way to revert my environment to using the in-built version of vpcSim() from the nlmixr2est package? Playing with namespaces directly is not something I know much about!
Thanks, O.
mattfidler
commented
4 months ago You would have to copy the old function to restore it. If you overwrote the function you would use:
vpcSim <- function(object, ..., keep=NULL, n=300,
pred=FALSE, seed=1009, nretry=50, minN=10,
normRelated=TRUE) {
checkmate::assertIntegerish(minN, len=1, any.missing=FALSE, lower=2)
checkmate::assertLogical(pred, len=1, any.missing=FALSE)
checkmate::assertIntegerish(nretry, len=1, any.missing=FALSE, lower=0)
checkmate::assertLogical(normRelated, len=1, any.missing=FALSE)
checkmate::assertCharacter(keep, null.ok=TRUE, pattern="^[.]*[a-zA-Z]+[a-zA-Z0-9._]*$")
checkmate::assertIntegerish(seed)
assignInMyNamespace(".finalUiCompressed", FALSE)
on.exit(assignInMyNamespace(".finalUiCompressed", TRUE))
set.seed(seed)
.si <- object$simInfo
.si$object <- eval(.getSimModel(object, hideIpred=FALSE))
.w <- which(names(.si) == "rx")
.si <- .si[-.w]
.si$nsim <- n
.si <- c(.si, list(...))
.pt <- proc.time()
.si$keep <- unique(c(keep, "nlmixrRowNums"))
.data <- .si$events
.data$nlmixrRowNums <- seq_along(.data[, 1])
if (normRelated) {
.ui <- object$ui
.predDf <-.ui$predDf
if (all(.predDf$dist %in% c("norm", "dnorm","t", "cauchy"))) {
} else {
if (is.null(.data$CMT)) {
.ds <- object$dataSav
.data$nlmixrRowNums <- seq_along(.data[,1])
.ds <- .ds[, c("CMT", "nlmixrRowNums")]
.data <- merge(.data, .ds, by ="nlmixrRowNums")
.data <- .data[order(.data$nlmixrRowNums),]
}
.lst <- .Call(`_nlmixr2est_filterNormalLikeAndDoses`,
.data$CMT, .predDf$distribution, .predDf$cmt)
.lst$nlmixrRowNums <- .data[.lst$filter, "nlmixrRowNums"]
if (.lst$nnorm == 0L) {
stop("need normal data for vpcSim (or use normRelated=FALSE)")
}
.data <- .data[.lst$filter, ]
}
}
.si$events <- .data
.si$thetaMat <- NULL
.si$dfSub <- NULL
.si$dfObs <- NULL
.si$returnType <- "data.frame.TBS"
.sim <- do.call(rxode2::rxSolve, .si)
# now look for how many have missing values
.w <- which(is.na(.sim$ipred))
.nretry <- 0
while (length(.w) > 0 && .nretry < nretry) {
.w <- which(is.na(.sim$ipred))
.simIds <- unique(.sim$sim.id[.w])
.sim <- .sim[!(.sim$sim.id %in% .simIds), ]
.sim$sim.id <- as.integer(factor(.sim$sim.id))
.mx <- max(.sim$sim.id)
.si$nsim <- n - .mx
.adjust <- FALSE
if (.si$nsim < minN) {
.si$nsim <- minN
.adjust <- TRUE
}
.sim2 <- do.call(rxode2::rxSolve, .si)
if (.adjust) {
# Select simulations without NA ipreds in them
.w <- which(is.na(.sim2$ipred))
.simIds <- unique(.sim2$sim.id[.w])
.allSimIds <- unique(.sim2$sim.id)
.simIds <- .allSimIds[!(.allSimIds %in% .simIds)]
.simIds <- .simIds[seq_len(min(n - .mx, length(.simIds)))]
.sim2 <- .sim2[.sim2$sim.id %in% .simIds, ]
}
.sim2$sim.id <- .sim2$sim.id + .mx
.sim <- rbind(.sim, .sim2)
.w <- which(is.na(.sim$ipred))
.nretry <- .nretry + 1
}
if (.nretry != 0) {
if (length(.w) == 0) {
warning("'NA' values in vpc or npde simulation, re-simulated until all simulations were successful",
call.=FALSE)
} else {
warning("'NA' values in vpc or npde simulation",
call.=FALSE)
}
}
if (pred) {
.si$nsim <- n # restore for pred
.si2 <- .si
.si2$params <- c(
.si$params, setNames(rep(0, dim(.si$omega)[1]), dimnames(.si$omega)[[2]]),
setNames(rep(0, dim(.si$sigma)[1]), dimnames(.si$sigma)[[2]]))
.si2$omega <- NULL
.si2$sigma <- NULL
.si2$returnType <- "data.frame"
.si2$nStud <- 1
.si2$nsim <- NULL
assignInMyNamespace(".lastPredSimulationInfo", .si2)
.sim2 <- do.call(rxode2::rxSolve, .si2)
.sim$pred <- .sim2$sim
}
.sim <- vpcNameDataCmts(object, .sim)
.cls <- c("nlmixr2vpcSim", class(.sim))
.fit <- object
.cls0 <- c("rxHidden", class(.fit))
attr(.cls0, ".foceiEnv") <- attr(class(.fit), ".foceiEnv")
class(.fit) <- .cls0
attr(.cls, "fit") <- .fit
class(.sim) <- .cls
return(.sim)
}
environment(vpcSim) <- asNamespace('nlmixr2est')
assignInNamespace("vpcSim", vpcSim, ns = "nlmixr2est")This is at the time of this issue submission.
However, it seems that there is an uncorrectable error with the simulation as it stands.
Here is the situation:
NA which implies that the simulation was performed improperly.Once there are NAs in the simulation, the problem cannot continue. This could be because of the sampling.
I cannot go any further without the data and model (sorry to say).
I do have some suggestions to get over this by trial and error:
addTable(vpc_tshoot_2, table = tableControl(cwres = TRUE, npde = TRUE, seed=472))Changing this seed may get over the problem you are having, but it is not guaranteed that this will fix the problem for you.
DrOJG
commented
4 months ago Hi Matt, Thank you for your help with this - a shame we couldn't get to the bottom of it for the moment but I can at least continue work now by either removing the npde from the table or removing the offending covariate term.
I will discuss with colleagues about potentially sharing more info, if not we may be able to revisit the discussion at a later date once we submit.
Enjoy the rest of your week! O.
mattfidler
commented
4 months ago This likely due to the platform independent random number generator implemented. This makes mac and windows share the same random numbers.
Unfortunately, mac and windows will not have the same random numbers with the new version as with the prior version.
DrOJG
commented
4 months ago I'm actually running the latest Ubuntu LTS version (odd the top of my head 22.04) if it makes a difference
mattfidler
commented
2 months ago I'm going to close this issue for now. If you can provide a reproducible example, I'd like to take another look at this issue and resolve it.
Hello there,
A piece of model code that was previously working for estimation with nlmixr2() is now no longer returning the data frame portion of the model fit object.
While estimation still works, and other portions of the model fit are recoverable (OFV, parameter estimates etc.), the main data frame is not there. Under run info, I see it now states "error calculating tables, returning without table step".
This seemed to coincide with some system updates, so I am wondering if there are any incompatibilities with new versions of compilers or other dependencies I should be aware of?
Many thanks, O.