Uncertainty ranges for low, md, and high values are not correctly created (at least in my example) · Issue #35 · certara/tidyvpc (original) (raw)

In the example below, all simulated values are reasonably good, but the simulated results are far from the expected values. My guess is that there is an error creating the bins in the simulated data.

suppressPackageStartupMessages(library(tidyvpc, quietly = TRUE)) library(nlmixr2, quietly = TRUE)

one_compartment <- function() { ini({ tka <- log(1.57); label("Ka") tcl <- log(2.72); label("Cl") tv <- log(31.5); label("V") eta_ka ~ 0.6 eta_cl ~ 0.3 eta_v ~ 0.1 add_sd <- 0.7 })

and a model block with the error specification and model specification

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) }) }

data_model <- theo_sd data_model$WTSTRATA <- ifelse(data_model$WT < median(data_model$WT), "Low weight", "High weight")

fit <- nlmixr2(one_compartment, data_model, est="saem", saemControl(print=0)) #> → loading into symengine environment... #> → pruning branches (if/else) of saem model... #> ✔ done #> → finding duplicate expressions in saem model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> → optimizing duplicate expressions in saem model... #> [====|====|====|====|====|====|====|====|====|====] 0:00:00 #> ✔ done #> rxode2 2.0.13.9000 using 8 threads (see ?getRxThreads) #> no cache: create with rxCreateCache() #> Calculating covariance matrix #> → loading into symengine environment... #> → pruning branches (if/else) of saem model... #> ✔ done #> → finding duplicate expressions in saem predOnly model 0... #> → finding duplicate expressions in saem predOnly model 1... #> → optimizing duplicate expressions in saem predOnly model 1... #> → finding duplicate expressions in saem predOnly model 2... #> ✔ done #> → Calculating residuals/tables #> ✔ done #> → compress origData in nlmixr2 object, save 7288 #> → compress phiM in nlmixr2 object, save 64048 #> → compress parHist in nlmixr2 object, save 9760 #> → compress saem0 in nlmixr2 object, save 30704

fit_vpcsim <- vpcSim(object = fit, keep = "WTSTRATA")

vpc <- observed(data_model, x=TIME, y=DV) %>% simulated(fit_vpcsim, x=time, y=sim) %>% stratify(~ WTSTRATA) %>% binning(bin = "jenks") %>% vpcstats()

vpc$stats #> WTSTRATA bin xbin qname y lo md hi #> 1: High weight [0,0.63) 0.000 q0.05 0.0000 -0.5566503 0.2496746 1.331612 #> 2: High weight [0,0.63) 0.000 q0.5 0.7400 3.8921545 5.3157246 6.796549 #> 3: High weight [0,0.63) 0.000 q0.95 7.2290 7.4868713 9.0411919 10.989930 #> 4: High weight [0.63,2.13) 1.135 q0.05 6.3265 -0.6620055 0.4778264 2.732001 #> 5: High weight [0.63,2.13) 1.135 q0.5 8.0000 3.9236263 5.5316531 6.754872 #> 6: High weight [0.63,2.13) 1.135 q0.95 9.9540 7.0465457 8.7930168 11.017270 #> 7: High weight [2.13,3.82) 3.530 q0.05 5.5690 -0.6199534 1.0852235 4.470570 #> 8: High weight [2.13,3.82) 3.530 q0.5 6.8500 2.6005538 5.2621121 7.691324 #> 9: High weight [2.13,3.82) 3.530 q0.95 8.1280 5.8337147 8.2528641 10.542600 #> 10: High weight [3.82,5.1) 5.020 q0.05 5.1590 -0.3739470 0.9844942 2.876889 #> 11: High weight [3.82,5.1) 5.020 q0.5 6.0800 3.6105409 5.2499374 6.915303 #> 12: High weight [3.82,5.1) 5.020 q0.95 8.0700 6.4690981 8.1877277 10.308895 #> 13: High weight [5.1,7.17) 7.030 q0.05 4.2330 -0.6383082 0.9801091 3.565333 #> 14: High weight [5.1,7.17) 7.030 q0.5 5.4000 2.6738969 5.2907440 7.170825 #> 15: High weight [5.1,7.17) 7.030 q0.95 8.0930 5.7688873 7.8856561 10.509108 #> 16: High weight [7.17,9.22) 9.000 q0.05 4.1490 -0.6912736 0.9872804 4.184987 #> 17: High weight [7.17,9.22) 9.000 q0.5 4.5700 3.3223763 5.2311310 7.037905 #> 18: High weight [7.17,9.22) 9.000 q0.95 6.4220 5.9220748 8.0771765 10.445583 #> 19: High weight [9.22,12.2) 12.000 q0.05 2.8460 -0.5199907 0.9724414 3.633877 #> 20: High weight [9.22,12.2) 12.000 q0.5 3.1600 3.3775260 5.2110285 7.159039 #> 21: High weight [9.22,12.2) 12.000 q0.95 5.4150 6.1720697 8.1718801 10.498052 #> 22: High weight [12.2,24.6] 24.235 q0.05 0.9070 -0.5749832 0.9416357 3.906557 #> 23: High weight [12.2,24.6] 24.235 q0.5 1.1350 2.5640919 5.1491458 6.968943 #> 24: High weight [12.2,24.6] 24.235 q0.95 3.5530 6.5890284 8.4005849 10.502630 #> 25: Low weight [0,1.02) 0.250 q0.05 0.0000 -0.7821223 0.2710994 1.493815 #> 26: Low weight [0,1.02) 0.250 q0.5 1.2500 3.9637494 5.3436445 6.676958 #> 27: Low weight [0,1.02) 0.250 q0.95 7.9820 7.2208578 8.8325933 10.707706 #> 28: Low weight [1.02,2.05) 1.990 q0.05 5.3675 -0.4897079 1.1370327 4.473624 #> 29: Low weight [1.02,2.05) 1.990 q0.5 6.6950 3.2127944 5.3481969 7.463388 #> 30: Low weight [1.02,2.05) 1.990 q0.95 9.6225 5.7939387 8.2028676 10.570098 #> 31: Low weight [2.05,5.07) 3.575 q0.05 5.5125 -0.6938220 0.6628048 3.911220 #> 32: Low weight [2.05,5.07) 3.575 q0.5 7.6950 3.1949944 5.2094577 7.385297 #> 33: Low weight [2.05,5.07) 3.575 q0.95 10.0030 6.0166708 8.2223557 10.670607 #> 34: Low weight [5.07,7.08) 7.020 q0.05 4.6100 -0.6466089 1.7813293 4.953796 #> 35: Low weight [5.07,7.08) 7.020 q0.5 6.5900 1.1325489 5.3456813 7.855727 #> 36: Low weight [5.07,7.08) 7.020 q0.95 8.2740 5.5283419 7.8334008 10.168535 #> 37: Low weight [7.08,9.38) 9.030 q0.05 3.7740 -0.4946345 1.3240811 4.811784 #> 38: Low weight [7.08,9.38) 9.030 q0.5 5.9000 2.6530119 5.1852825 7.431180 #> 39: Low weight [7.08,9.38) 9.030 q0.95 7.6380 5.7928305 8.0354830 10.563015 #> 40: Low weight [9.38,12.1) 12.050 q0.05 3.6980 -0.4840723 1.3381509 4.455000 #> 41: Low weight [9.38,12.1) 12.050 q0.5 4.5700 2.4324854 5.3852527 7.784143 #> 42: Low weight [9.38,12.1) 12.050 q0.95 6.8480 5.5211620 8.1486792 11.071007 #> 43: Low weight [12.1,24.4] 24.115 q0.05 0.9325 -0.7706144 1.4873224 4.483279 #> 44: Low weight [12.1,24.4] 24.115 q0.5 1.3700 2.1159037 5.1524565 7.671587 #> 45: Low weight [12.1,24.4] 24.115 q0.95 2.6225 5.8568803 8.1987823 10.677722 #> WTSTRATA bin xbin qname y lo md hi plot(vpc)