mattfidler
commented
5 months ago This should work is the development version of nlmixr2est. I'm assuming you are not using it. There is a note about how fixing was changed in that version.
If you are using the development nlmixr2 let me know and I will look into this further. Otherwise I will continue working toward releasing this on CRAN.
If you want an explanation of the differences I can also provide one for you.
Best,
Matt
On Mon, May 20, 2024, 10:35 PM Zhonghui Huang @.***> wrote:
Hi @mattfidler https://github.com/mattfidler,
I tried to introduce covariate models with fixed power values by fix() function in the ini block, however, I find that the fix() function might not work if we use saem not focei, I attached my test results here by using the theo_sd dataset.
Test 1. Covariate model format: ( covwt<- fix(0.01) and use covwt lwt, where lwt = log(theo_sd2$WT/70)) by saem Test 2. Covariate model format: ( covwt<- fix(0.75) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by saem Test 3. Covariate model format: ( covwt<- fix(100) and use covwt*lwt, where lwt = log(theo_sd2$WT/70)) by saem
Test 4. Covariate model format: ( covwt<- fix(0.01) and use covwt lwt, where lwt = log(theo_sd2$WT/70)) by focei Test 5. Covariate model format: ( covwt<- fix(0.75) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by focei Test 6. Covariate model format: ( covwt<- fix(100) and use covwt*lwt, where lwt = log(theo_sd2$WT/70)) by focei
Test1-3 shared totally same parameter estimates in each iteration, and share the same final parameter estimate results. But they shouldnt be same given the different setting of power, especially the first and third are extreme values.
However, apparently, focei can work. When focei was used in test 4, test 5 and test 6, we can see the change of parameter estimates between iterations and the model estimate results are also different.
So, if we would like to introduce a fixed allometric scalling with a power of 0.75 and 1 and also choose saem, currently, my option is to write the equation directly in the model block, e.g.,
cl <- exp(tcl + eta.cl http://eta.cl + 0.75 log(WT/70)) cl<-exp(tcl + eta.cl http://eta.cl) (WT/70) ^0.75**
(BTW, cl=exp(tcl+eta.cl+lwt*0.75),lwt = log(theo_sd2$WT/70)), also doesnt work for me) But still hope the fix function can work then can see the fixed values in the nlmixr2 output.
Test1. Use fix() function to introduce a fixed covariate model, extreme power 0.01
Use fix() function to introduce a fixed covariate model> #==============================================================================#> > library(nlmixr2)> theo_sd2<-theo_sd> theo_sd2$lwt<-log(theo_sd2$WT/70)> > ## The basic model consists of an ini block that has initial estimates> one.compartment <- function() {+ ini({+ tka <- log(1.57); label("Ka")+ tcl <- log(2.72); label("Cl")+ tv <- log(31.5); label("V")+ covwt<- fix(0.01);+ + 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 + covwtlwt)+ 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)+ })+ }> > ## The fit is performed by the function nlmixr/nlmixr2 specifying the model, data and estimate> fit0<- nlmixr2(one.compartment, theo_sd2, est="saem", saemControl(print=10,seed = 1234))
→ 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
✔ doneparams: tka tcl tv V(eta.ka) V(eta.cl) V(eta.v) add.sd001: 0.426017 1.083247 3.539728 0.592610 0.285000 0.095000 2.408191 010: 0.501829 0.967554 3.476164 0.410733 0.179621 0.059874 0.731773 020: 0.454941 1.014119 3.459420 0.346308 0.107546 0.035849 0.705265 030: 0.477338 0.995463 3.474201 0.394533 0.069945 0.021464 0.704248 040: 0.439342 1.028859 3.449393 0.431745 0.085687 0.016491 0.681260 050: 0.380047 1.066219 3.426244 0.415438 0.077686 0.010803 0.707764 060: 0.412688 1.004636 3.438032 0.424521 0.091897 0.012759 0.706832 070: 0.452603 1.009378 3.463855 0.360588 0.075677 0.017905 0.703528 080: 0.370318 1.046738 3.420545 0.397710 0.081115 0.015644 0.711544 090: 0.460203 1.023362 3.449975 0.440390 0.071412 0.023928 0.693081 100: 0.450080 1.018632 3.455627 0.415294 0.067890 0.022319 0.687304 110: 0.455493 1.027576 3.441667 0.445716 0.075984 0.016721 0.674659 120: 0.468103 1.016301 3.451920 0.452389 0.075484 0.018345 0.685905 130: 0.471006 0.991160 3.469804 0.360780 0.074823 0.023260 0.711663 140: 0.437146 1.012946 3.446918 0.418566 0.063419 0.029043 0.685813 150: 0.486432 0.992229 3.473682 0.399591 0.063687 0.027025 0.697351 160: 0.482319 1.013383 3.462588 0.413951 0.092843 0.021622 0.678935 170: 0.503791 1.027372 3.468940 0.483595 0.058152 0.020509 0.689687 180: 0.439974 1.007908 3.454167 0.371606 0.074960 0.014613 0.678929 190: 0.427293 1.037723 3.439008 0.374604 0.062022 0.013980 0.692228 200: 0.440119 1.040741 3.458134 0.355890 0.074055 0.013313 0.721150 210: 0.452073 1.020168 3.446681 0.385688 0.072180 0.015044 0.700821 220: 0.447489 1.020724 3.448970 0.383226 0.073030 0.015859 0.699757 230: 0.444974 1.020360 3.447842 0.384076 0.073989 0.015420 0.700211 240: 0.445238 1.022240 3.447344 0.383898 0.074468 0.015518 0.702492 250: 0.442763 1.023440 3.447074 0.384727 0.073796 0.015493 0.700587 260: 0.444879 1.022585 3.447789 0.388161 0.072947 0.015804 0.699259 270: 0.442805 1.021528 3.446868 0.393174 0.072693 0.016076 0.697323 280: 0.447570 1.018403 3.448971 0.394812 0.073098 0.015982 0.697241 290: 0.451412 1.015877 3.450195 0.391768 0.072871 0.016015 0.697483 300: 0.449270 1.015384 3.450110 0.391312 0.071840 0.016510 0.697923 310: 0.447051 1.015763 3.449151 0.389408 0.072454 0.016442 0.697973 320: 0.449876 1.014353 3.449954 0.391474 0.072775 0.016388 0.698833 330: 0.451111 1.013431 3.450882 0.391544 0.072851 0.016440 0.698968 340: 0.450308 1.014217 3.450329 0.392854 0.072530 0.016576 0.699140 350: 0.448806 1.014773 3.449905 0.392169 0.072751 0.016395 0.699271 360: 0.449157 1.015164 3.449856 0.391511 0.072889 0.016301 0.699652 370: 0.447701 1.015830 3.449027 0.393459 0.072397 0.016491 0.699513 380: 0.448015 1.015419 3.449165 0.394348 0.072100 0.016563 0.699113 390: 0.447556 1.015096 3.449327 0.394829 0.071962 0.016632 0.698792 400: 0.448474 1.014165 3.449721 0.395529 0.071723 0.016774 0.698800 410: 0.448255 1.014056 3.449827 0.393472 0.071843 0.016773 0.699147 420: 0.448896 1.014260 3.450057 0.393950 0.072037 0.016694 0.699327 430: 0.449246 1.014099 3.450090 0.393667 0.071786 0.016717 0.699200 440: 0.449055 1.013679 3.450186 0.391934 0.071733 0.016784 0.699169 450: 0.449796 1.013347 3.450530 0.393240 0.071856 0.016882 0.699053 460: 0.450259 1.013232 3.450704 0.393598 0.071889 0.016968 0.698976 470: 0.450742 1.012861 3.451049 0.393049 0.072178 0.016902 0.698485 480: 0.451073 1.012971 3.451105 0.393940 0.071809 0.016975 0.698098 490: 0.450468 1.013045 3.451226 0.393096 0.071582 0.017017 0.698059 500: 0.450520 1.013017 3.451320 0.392231 0.071404 0.017024 0.698030 Calculating covariance matrix [====|====|====|====|====|====|====|====|====|====] 0:00:00
→ loading into symengine environment... → pruning branches (
if/else) of saem model... ✔ done → finding duplicate expressions in saem predOnly model 0... [====|====|====|====|====|====|====|====|====|====] 0:00:00→ finding duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00
→ optimizing duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00
→ finding duplicate expressions in saem predOnly model 2... [====|====|====|====|====|====|====|====|====|====] 0:00:00
✔ done → Calculating residuals/tables ✔ done → compress origData in nlmixr2 object, save 7056 → compress phiM in nlmixr2 object, save 82136 → compress parHistData in nlmixr2 object, save 13888 → compress saem0 in nlmixr2 object, save 30864> print(fit0) ── nlmixr² SAEM OBJF by FOCEi approximation ──
Gaussian/Laplacian Likelihoods: AIC(fit0) or fit0$objf etc. FOCEi CWRES & Likelihoods: addCwres(fit0)
── Time (sec fit0$time): ──
setup covariance saem table compress otherelapsed 0.000658 0.005003 1.944 0.034 0.017 0.581339── Population Parameters (fit0$parFixed or fit0$parFixedDf): ──
Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%tka Ka 0.451 0.191 42.4 1.57 (1.08, 2.28) 69.3 -0.609% tcl Cl 1.01 0.0844 8.33 2.75 (2.33, 3.25) 27.2 4.13% tv V 3.45 0.0448 1.3 31.5 (28.9, 34.4) 13.1 9.48% covwt 0.01 FIXED FIXED 0.01 add.sd 0.698 0.698Covariance Type (fit0$covMethod): linFim No correlations in between subject variability (BSV) matrix Full BSV covariance (fit0$omega) or correlation (fit0$omegaR; diagonals=SDs) Distribution stats (mean/skewness/kurtosis/p-value) available in fit0$shrink Censoring (fit0$censInformation): No censoring
── Fit Data (object fit0 is a modified tibble): ──# A tibble: 132 × 20 ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp depot center ka cl v tad dosenum lwt
1 1 0 0.74 0 0.74 0 0.74 1.06 0.104 -0.482 -0.0740 0 320. 0 1.74 1.70 29.3 0 1 0.1292 1 0.25 2.84 3.25 -0.414 3.83 -0.987 -1.41 0.104 -0.482 -0.0740 3.83 207. 112. 1.74 1.70 29.3 0.25 1 0.1293 1 0.57 6.57 5.83 0.741 6.75 -0.176 -0.252 0.104 -0.482 -0.0740 6.75 119. 198. 1.74 1.70 29.3 0.57 1 0.129# ℹ 129 more rows# ℹ Use `print(n = ...)` to see more rows *Test2. Use fix() function to introduce a fixed covariate model, 0.75* > > library(nlmixr2)> theo_sd2<-theo_sd> theo_sd2$lwt<-log(theo_sd2$WT/70)> > ## The basic model consists of an ini block that has initial estimates> one.compartment <- function() {+ ini({+ tka <- log(1.57); label("Ka")+ tcl <- log(2.72); label("Cl")+ tv <- log(31.5); label("V")+ covwt<- fix(0.75);+ + 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 + covwt*lwt)+ 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)+ })+ }> > ## The fit is performed by the function nlmixr/nlmixr2 specifying the model, data and estimate> fit <- nlmixr2(one.compartment, theo_sd2, est="saem", saemControl(print=10,seed = 1234)) → 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 ✔ doneparams: tka tcl tv V(eta.ka) V(eta.cl) V(eta.v) add.sd001: 0.426017 1.083247 3.539728 0.592610 0.285000 0.095000 2.408191 010: 0.501829 0.967554 3.476164 0.410733 0.179621 0.059874 0.731773 020: 0.454941 1.014119 3.459420 0.346308 0.107546 0.035849 0.705265 030: 0.477338 0.995463 3.474201 0.394533 0.069945 0.021464 0.704248 040: 0.439342 1.028859 3.449393 0.431745 0.085687 0.016491 0.681260 050: 0.380047 1.066219 3.426244 0.415438 0.077686 0.010803 0.707764 060: 0.412688 1.004636 3.438032 0.424521 0.091897 0.012759 0.706832 070: 0.452603 1.009378 3.463855 0.360588 0.075677 0.017905 0.703528 080: 0.370318 1.046738 3.420545 0.397710 0.081115 0.015644 0.711544 090: 0.460203 1.023362 3.449975 0.440390 0.071412 0.023928 0.693081 100: 0.450080 1.018632 3.455627 0.415294 0.067890 0.022319 0.687304 110: 0.455493 1.027576 3.441667 0.445716 0.075984 0.016721 0.674659 120: 0.468103 1.016301 3.451920 0.452389 0.075484 0.018345 0.685905 130: 0.471006 0.991160 3.469804 0.360780 0.074823 0.023260 0.711663 140: 0.437146 1.012946 3.446918 0.418566 0.063419 0.029043 0.685813 150: 0.486432 0.992229 3.473682 0.399591 0.063687 0.027025 0.697351 160: 0.482319 1.013383 3.462588 0.413951 0.092843 0.021622 0.678935 170: 0.503791 1.027372 3.468940 0.483595 0.058152 0.020509 0.689687 180: 0.439974 1.007908 3.454167 0.371606 0.074960 0.014613 0.678929 190: 0.427293 1.037723 3.439008 0.374604 0.062022 0.013980 0.692228 200: 0.440119 1.040741 3.458134 0.355890 0.074055 0.013313 0.721150 210: 0.452073 1.020168 3.446681 0.385688 0.072180 0.015044 0.700821 220: 0.447489 1.020724 3.448970 0.383226 0.073030 0.015859 0.699757 230: 0.444974 1.020360 3.447842 0.384076 0.073989 0.015420 0.700211 240: 0.445238 1.022240 3.447344 0.383898 0.074468 0.015518 0.702492 250: 0.442763 1.023440 3.447074 0.384727 0.073796 0.015493 0.700587 260: 0.444879 1.022585 3.447789 0.388161 0.072947 0.015804 0.699259 270: 0.442805 1.021528 3.446868 0.393174 0.072693 0.016076 0.697323 280: 0.447570 1.018403 3.448971 0.394812 0.073098 0.015982 0.697241 290: 0.451412 1.015877 3.450195 0.391768 0.072871 0.016015 0.697483 300: 0.449270 1.015384 3.450110 0.391312 0.071840 0.016510 0.697923 310: 0.447051 1.015763 3.449151 0.389408 0.072454 0.016442 0.697973 320: 0.449876 1.014353 3.449954 0.391474 0.072775 0.016388 0.698833 330: 0.451111 1.013431 3.450882 0.391544 0.072851 0.016440 0.698968 340: 0.450308 1.014217 3.450329 0.392854 0.072530 0.016576 0.699140 350: 0.448806 1.014773 3.449905 0.392169 0.072751 0.016395 0.699271 360: 0.449157 1.015164 3.449856 0.391511 0.072889 0.016301 0.699652 370: 0.447701 1.015830 3.449027 0.393459 0.072397 0.016491 0.699513 380: 0.448015 1.015419 3.449165 0.394348 0.072100 0.016563 0.699113 390: 0.447556 1.015096 3.449327 0.394829 0.071962 0.016632 0.698792 400: 0.448474 1.014165 3.449721 0.395529 0.071723 0.016774 0.698800 410: 0.448255 1.014056 3.449827 0.393472 0.071843 0.016773 0.699147 420: 0.448896 1.014260 3.450057 0.393950 0.072037 0.016694 0.699327 430: 0.449246 1.014099 3.450090 0.393667 0.071786 0.016717 0.699200 440: 0.449055 1.013679 3.450186 0.391934 0.071733 0.016784 0.699169 450: 0.449796 1.013347 3.450530 0.393240 0.071856 0.016882 0.699053 460: 0.450259 1.013232 3.450704 0.393598 0.071889 0.016968 0.698976 470: 0.450742 1.012861 3.451049 0.393049 0.072178 0.016902 0.698485 480: 0.451073 1.012971 3.451105 0.393940 0.071809 0.016975 0.698098 490: 0.450468 1.013045 3.451226 0.393096 0.071582 0.017017 0.698059 500: 0.450520 1.013017 3.451320 0.392231 0.071404 0.017024 0.698030 Calculating covariance matrix [====|====|====|====|====|====|====|====|====|====] 0:00:00 → loading into symengine environment... → pruning branches (`if`/`else`) of saem model... ✔ done → finding duplicate expressions in saem predOnly model 0... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → finding duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → optimizing duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → finding duplicate expressions in saem predOnly model 2... [====|====|====|====|====|====|====|====|====|====] 0:00:00 ✔ done → Calculating residuals/tables ✔ done → compress origData in nlmixr2 object, save 7056 → compress phiM in nlmixr2 object, save 94688 → compress parHistData in nlmixr2 object, save 13888 → compress saem0 in nlmixr2 object, save 30896> print(fit) ── nlmixr² SAEM OBJF by FOCEi approximation ── Gaussian/Laplacian Likelihoods: AIC(fit) or fit$objf etc. FOCEi CWRES & Likelihoods: addCwres(fit) ── Time (sec fit$time): ── setup covariance saem table compress otherelapsed 0.000618 0.003003 1.911 0.035 0.018 0.510379 ── Population Parameters (fit$parFixed or fit$parFixedDf): ── Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%tka Ka 0.451 0.191 42.4 1.57 (1.08, 2.28) 69.3 -0.609% tcl Cl 1.01 0.0844 8.33 2.75 (2.33, 3.25) 27.2 4.13% tv V 3.45 0.0448 1.3 31.5 (28.9, 34.4) 13.1 9.48% covwt 0.75 FIXED FIXED 0.75 add.sd 0.698 0.698 Covariance Type (fit$covMethod): linFim No correlations in between subject variability (BSV) matrix Full BSV covariance (fit$omega) or correlation (fit$omegaR; diagonals=SDs) Distribution stats (mean/skewness/kurtosis/p-value) available in fit$shrink Censoring (fit$censInformation): No censoring ── Fit Data (object fit is a modified tibble): ──# A tibble: 132 × 20 ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp depot center ka cl v tad dosenum lwt 1 1 0 0.74 0 0.74 0 0.74 1.06 0.104 -0.482 -0.0740 0 320. 0 1.74 1.87 29.3 0 1 0.1292 1 0.25 2.84 3.25 -0.410 3.82 -0.984 -1.41 0.104 -0.482 -0.0740 3.82 207. 112. 1.74 1.87 29.3 0.25 1 0.1293 1 0.57 6.57 5.81 0.758 6.73 -0.163 -0.234 0.104 -0.482 -0.0740 6.73 119. 197. 1.74 1.87 29.3 0.57 1 0.129 *Test 3. Use fix() function to introduce a fixed covariate model, use an extreme value, 100* > library(nlmixr2)> theo_sd2<-theo_sd> theo_sd2$lwt<-log(theo_sd2$WT/70)> > > ## The basic model consists of an ini block that has initial estimates> one.compartment <- function() {+ ini({+ tka <- log(1.57); label("Ka")+ tcl <- log(2.72); label("Cl")+ tv <- log(31.5); label("V")+ covwt<- fix(100);+ + 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 + covwt*lwt)+ 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)+ })+ }> > ## The fit is performed by the function nlmixr/nlmixr2 specifying the model, data and estimate> fit2<- nlmixr2(one.compartment, theo_sd2, est="saem", saemControl(print=10,seed = 1234)) → 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 ✔ doneparams: tka tcl tv V(eta.ka) V(eta.cl) V(eta.v) add.sd001: 0.426017 1.083247 3.539728 0.592610 0.285000 0.095000 2.408191 010: 0.501829 0.967554 3.476164 0.410733 0.179621 0.059874 0.731773 020: 0.454941 1.014119 3.459420 0.346308 0.107546 0.035849 0.705265 030: 0.477338 0.995463 3.474201 0.394533 0.069945 0.021464 0.704248 040: 0.439342 1.028859 3.449393 0.431745 0.085687 0.016491 0.681260 050: 0.380047 1.066219 3.426244 0.415438 0.077686 0.010803 0.707764 060: 0.412688 1.004636 3.438032 0.424521 0.091897 0.012759 0.706832 070: 0.452603 1.009378 3.463855 0.360588 0.075677 0.017905 0.703528 080: 0.370318 1.046738 3.420545 0.397710 0.081115 0.015644 0.711544 090: 0.460203 1.023362 3.449975 0.440390 0.071412 0.023928 0.693081 100: 0.450080 1.018632 3.455627 0.415294 0.067890 0.022319 0.687304 110: 0.455493 1.027576 3.441667 0.445716 0.075984 0.016721 0.674659 120: 0.468103 1.016301 3.451920 0.452389 0.075484 0.018345 0.685905 130: 0.471006 0.991160 3.469804 0.360780 0.074823 0.023260 0.711663 140: 0.437146 1.012946 3.446918 0.418566 0.063419 0.029043 0.685813 150: 0.486432 0.992229 3.473682 0.399591 0.063687 0.027025 0.697351 160: 0.482319 1.013383 3.462588 0.413951 0.092843 0.021622 0.678935 170: 0.503791 1.027372 3.468940 0.483595 0.058152 0.020509 0.689687 180: 0.439974 1.007908 3.454167 0.371606 0.074960 0.014613 0.678929 190: 0.427293 1.037723 3.439008 0.374604 0.062022 0.013980 0.692228 200: 0.440119 1.040741 3.458134 0.355890 0.074055 0.013313 0.721150 210: 0.452073 1.020168 3.446681 0.385688 0.072180 0.015044 0.700821 220: 0.447489 1.020724 3.448970 0.383226 0.073030 0.015859 0.699757 230: 0.444974 1.020360 3.447842 0.384076 0.073989 0.015420 0.700211 240: 0.445238 1.022240 3.447344 0.383898 0.074468 0.015518 0.702492 250: 0.442763 1.023440 3.447074 0.384727 0.073796 0.015493 0.700587 260: 0.444879 1.022585 3.447789 0.388161 0.072947 0.015804 0.699259 270: 0.442805 1.021528 3.446868 0.393174 0.072693 0.016076 0.697323 280: 0.447570 1.018403 3.448971 0.394812 0.073098 0.015982 0.697241 290: 0.451412 1.015877 3.450195 0.391768 0.072871 0.016015 0.697483 300: 0.449270 1.015384 3.450110 0.391312 0.071840 0.016510 0.697923 310: 0.447051 1.015763 3.449151 0.389408 0.072454 0.016442 0.697973 320: 0.449876 1.014353 3.449954 0.391474 0.072775 0.016388 0.698833 330: 0.451111 1.013431 3.450882 0.391544 0.072851 0.016440 0.698968 340: 0.450308 1.014217 3.450329 0.392854 0.072530 0.016576 0.699140 350: 0.448806 1.014773 3.449905 0.392169 0.072751 0.016395 0.699271 360: 0.449157 1.015164 3.449856 0.391511 0.072889 0.016301 0.699652 370: 0.447701 1.015830 3.449027 0.393459 0.072397 0.016491 0.699513 380: 0.448015 1.015419 3.449165 0.394348 0.072100 0.016563 0.699113 390: 0.447556 1.015096 3.449327 0.394829 0.071962 0.016632 0.698792 400: 0.448474 1.014165 3.449721 0.395529 0.071723 0.016774 0.698800 410: 0.448255 1.014056 3.449827 0.393472 0.071843 0.016773 0.699147 420: 0.448896 1.014260 3.450057 0.393950 0.072037 0.016694 0.699327 430: 0.449246 1.014099 3.450090 0.393667 0.071786 0.016717 0.699200 440: 0.449055 1.013679 3.450186 0.391934 0.071733 0.016784 0.699169 450: 0.449796 1.013347 3.450530 0.393240 0.071856 0.016882 0.699053 460: 0.450259 1.013232 3.450704 0.393598 0.071889 0.016968 0.698976 470: 0.450742 1.012861 3.451049 0.393049 0.072178 0.016902 0.698485 480: 0.451073 1.012971 3.451105 0.393940 0.071809 0.016975 0.698098 490: 0.450468 1.013045 3.451226 0.393096 0.071582 0.017017 0.698059 500: 0.450520 1.013017 3.451320 0.392231 0.071404 0.017024 0.698030 Calculating covariance matrix [====|====|====|====|====|====|====|====|====|====] 0:00:00 → loading into symengine environment... → pruning branches (`if`/`else`) of saem model... ✔ done → finding duplicate expressions in saem predOnly model 0... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → finding duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → optimizing duplicate expressions in saem predOnly model 1... [====|====|====|====|====|====|====|====|====|====] 0:00:00 → finding duplicate expressions in saem predOnly model 2... [====|====|====|====|====|====|====|====|====|====] 0:00:00 ✔ done → Calculating residuals/tables ✔ done → compress origData in nlmixr2 object, save 7056 → compress phiM in nlmixr2 object, save 139544 → compress parHistData in nlmixr2 object, save 13888 → compress saem0 in nlmixr2 object, save 30928> print(fit2) ── nlmixr² SAEM OBJF by FOCEi approximation ── Gaussian/Laplacian Likelihoods: AIC(fit2) or fit2$objf etc. FOCEi CWRES & Likelihoods: addCwres(fit2) ── Time (sec fit2$time): ── setup covariance saem table compress otherelapsed 0.000639 0.005003 1.893 0.032 0.016 0.519358 ── Population Parameters (fit2$parFixed or fit2$parFixedDf): ── Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%tka Ka 0.451 0.191 42.4 1.57 (1.08, 2.28) 69.3 -0.609% tcl Cl 1.01 0.0844 8.33 2.75 (2.33, 3.25) 27.2 4.13% tv V 3.45 0.0448 1.3 31.5 (28.9, 34.4) 13.1 9.48% covwt 100 FIXED FIXED 100 add.sd 0.698 0.698 Covariance Type (fit2$covMethod): linFim No correlations in between subject variability (BSV) matrix Full BSV covariance (fit2$omega) or correlation (fit2$omegaR; diagonals=SDs) Distribution stats (mean/skewness/kurtosis/p-value) available in fit2$shrink Censoring (fit2$censInformation): No censoring ── Fit Data (object fit2 is a modified tibble): ──# A tibble: 132 × 20 ID TIME DV PRED RES IPRED IRES IWRES eta.ka eta.cl eta.v cp depot center ka cl v tad dosenum lwt 1 1 0 0.74 0 0.74 0 0.74 1.06 0.104 -0.482 -0.0740 0 320. 0 1.74 648497. 29.3 0 1 0.1292 1 0.25 2.84 0.000323 2.84 0.000556 2.84 4.07 0.104 -0.482 -0.0740 0.000556 207. 0.0163 1.74 648497. 29.3 0.25 1 0.1293 1 0.57 6.57 0.000195 6.57 0.000318 6.57 9.41 0.104 -0.482 -0.0740 0.000318 119. 0.00933 1.74 648497. 29.3 0.57 1 0.129# ℹ 129 more rows *Test4. Use fix() function to introduce a fixed covariate model, run by FOCEI* > library(nlmixr2)> theo_sd2<-theo_sd> theo_sd2$lwt<-log(theo_sd2$WT/70)> > ## The basic model consists of an ini block that has initial estimates> one.compartment <- function() {+ ini({+ tka <- log(1.57); label("Ka")+ tcl <- log(2.72); label("Cl")+ tv <- log(31.5); label("V")+ covwt<- fix(0.01);+ + 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 + covwt*lwt)+ 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)+ })+ }> > ## The fit is performed by the function nlmixr/nlmixr2 specifying the model, data and estimate> fit3 <- nlmixr2(one.compartment, theo_sd2, est="focei", foceiControl(print=10,seed = 1234))Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximationF: Forward difference gradient approximationC: Central difference gradient approximationM: Mixed forward and central difference gradient approximationUnscaled parameters for Omegas=chol(solve(omega));Diagonals are transformed, as specified by foceiControl(diagXform=)|-----+---------------+-----------+-----------+-----------+-----------+-----------+-----------+-----------|| #| Objective Fun | tka | tcl | tv | add.sd | o1 | o2 | o3 || 10| 122.02628 | -0.9784 | -0.6251 | 0.9768 | -0.8380 | -0.3995 | 0.1663 | 0.3299 || U| 122.02628 | 0.4726 | 1.009 | 3.427 | 0.6943 | 1.186 | 1.770 | 2.028 || X| 122.02628 | 1.604 | 2.743 | 30.78 | 0.6943 | 1.186 | 1.770 | 2.028 || F| Forward Diff. | -0.3764 | 1.979 | 14.59 | 19.36 | -1.889 | -5.481 | -6.625 || 20| 117.40915 | -0.9623 | -0.6467 | 1.019 | -0.8353 | -0.2340 | 0.2912 | 1.147 || U| 117.40915 | 0.4888 | 0.9874 | 3.469 | 0.6982 | 1.244 | 1.863 | 2.488 || X| 117.40915 | 1.630 | 2.684 | 32.12 | 0.6982 | 1.244 | 1.863 | 2.488 || 30| 117.12000 | -0.9426 | -0.6200 | 1.016 | -0.8331 | -0.2607 | 0.3311 | 1.251 || U| 117.12 | 0.5085 | 1.014 | 3.466 | 0.7013 | 1.235 | 1.892 | 2.546 || X| 117.12 | 1.663 | 2.757 | 32.00 | 0.7013 | 1.235 | 1.892 | 2.546 || F| Forward Diff. | 0.08544 | -1.469 | -0.06992 | 12.16 | -1.051 | -1.599 | -1.055 || 40| 116.83087 | -0.9818 | -0.6218 | 1.008 | -0.8379 | -0.2491 | 0.3738 | 1.409 || U| 116.83087 | 0.4692 | 1.012 | 3.458 | 0.6944 | 1.239 | 1.924 | 2.635 || X| 116.83087 | 1.599 | 2.752 | 31.75 | 0.6944 | 1.239 | 1.924 | 2.635 || F| Forward Diff. | -5.067 | -0.9498 | -0.7960 | 0.3804 | -0.2611 | 0.07740 | -0.1536 || 50| 116.78891 | -0.9863 | -0.6232 | 1.011 | -0.8371 | -0.2257 | 0.3941 | 1.469 || U| 116.78891 | 0.4648 | 1.011 | 3.461 | 0.6956 | 1.247 | 1.939 | 2.669 || X| 116.78891 | 1.592 | 2.748 | 31.83 | 0.6956 | 1.247 | 1.939 | 2.669 |calculating covariance matrix [====|====|====|====|====|====|====|====|====|====] 0:00:00 done → Calculating residuals/tables ✔ done → compress origData in nlmixr2 object, save 7056 → compress parHistData in nlmixr2 object, save 6400> print(fit3) ── nlmixr² FOCEi (outer: nlminb) ── OBJF AIC BIC Log-likelihood Condition#(Cov) Condition#(Cor)FOCEi 116.7875 373.3872 393.5669 -179.6936 63.01769 4.138121 ── Time (sec fit3$time): ── setup optimize covariance table compress otherelapsed 0.000659 0.222847 0.222848 0.022 0.008 1.186646 ── Population Parameters (fit3$parFixed or fit3$parFixedDf): ── Parameter Est. SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%tka Ka 0.464 0.226 48.7 1.59 (1.02, 2.48) 71.5 2.83% tcl Cl 1.01 0.368 36.4 2.75 (1.34, 5.65) 27.1 5.05% tv V 3.46 0.0581 1.68 31.9 (28.4, 35.7) 14.1 11.1% covwt 0.01 FIXED FIXED 0.01 add.sd 0.695 0.695 Covariance Type (fit3$covMethod): r,s No correlations in between subject variability (BSV) matrix Full BSV covariance (fit3$omega) or correlation (fit3$omegaR; diagonals=SDs) Distribution stats (mean/skewness/kurtosis/p-value) available in fit3$shrink Information about run found (fit3$runInfo): • gradient problems with initial estimate and covariance; see $scaleInfo • ETAs were reset to zero during optimization; (Can control by foceiControl(resetEtaP=.)) • initial ETAs were nudged; (can control by foceiControl(etaNudge=., etaNudge2=)) Censoring (fit3$censInformation): No censoring Minimization message (fit3$message): relative convergence (4) ── Fit Data (object fit3 is a modified tibble): ──# A tibble: 132 × 23 ID TIME DV PRED RES WRES IPRED IRES IWRES CPRED CRES CWRES eta.ka eta.cl eta.v depot center ka cl v tad dosenum lwt 1 1 0 0.74 0 0.74 1.06 0 0.74 1.06 0 0.74 1.06 0.0865 -0.476 -0.0924 320. 0 1.73 1.71 29.0 0 1 0.1292 1 0.25 2.84 3.26 -0.419 -0.221 3.85 -1.01 -1.45 3.21 -0.374 -0.173 0.0865 -0.476 -0.0924 207. 112. 1.73 1.71 29.0 0.25 1 0.1293 1 0.57 6.57 5.82 0.748 0.297 6.79 -0.215 -0.310 5.77 0.805 0.288 0.0865 -0.476 -0.0924 119. 197. 1.73 1.71 29.0 0.57 1 0.129# ℹ 129 more rows# ℹ Use `print(n = ...)` to see more rows *Test5. Use fix() function to introduce a fixed covariate model, run by FOCEI* > library(nlmixr2)> theo_sd2<-theo_sd> theo_sd2$lwt<-log(theo_sd2$WT/70)> > ## The basic model consists of an ini block that has initial estimates> one.compartment <- function() {+ ini({+ tka <- log(1.57); label("Ka")+ tcl <- log(2.72); label("Cl")+ tv
ZhonghuiHuang
Hi @mattfidler,
I tried to introduce covariate models with fixed power values by fix() function in the ini block, however, I find that the fix() function might not work if we use saem not focei, I attached my test results here by using the theo_sd dataset.
Test 1. Covariate model format: ( covwt<- fix(0.01) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by saem Test 2. Covariate model format: ( covwt<- fix(0.75) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by saem Test 3. Covariate model format: ( covwt<- fix(100) and use covwt*lwt, where lwt = log(theo_sd2$WT/70)) by saem
Test 4. Covariate model format: ( covwt<- fix(0.01) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by focei Test 5. Covariate model format: ( covwt<- fix(0.75) and use covwtlwt, where lwt = log(theo_sd2$WT/70)) by focei Test 6. Covariate model format: ( covwt<- fix(100) and use covwt*lwt, where lwt = log(theo_sd2$WT/70)) by focei
Test1-3 shared totally same parameter estimates in each iteration, and share the same final parameter estimate results. But they shouldnt be same given the different setting of power, especially the first and third are extreme values.
However, apparently, focei can work. When focei was used in test 4, test 5 and test 6, we can see the change of parameter estimates between iterations and the model estimate results are also different.
So, if we would like to introduce a fixed allometric scalling with a power of 0.75 and 1 and also choose saem, currently, my option is to write the equation directly in the model block, e.g.,
cl <- exp(tcl + eta.cl + 0.75 log(WT/70)) cl<-exp(tcl + eta.cl) (WT/70) ^0.75
(BTW, cl=exp(tcl+eta.cl+lwt*0.75),lwt = log(theo_sd2$WT/70)), also doesnt work for me) But still hope the fix function can work then can see the fixed values in the nlmixr2 output.
Test1. Use fix() function to introduce a fixed covariate model, extreme power 0.01
Test2. Use fix() function to introduce a fixed covariate model, 0.75
Test 3. Use fix() function to introduce a fixed covariate model, use an extreme value, 100
Test4. Use fix() function to introduce a fixed covariate model, run by FOCEI
Test5. Use fix() function to introduce a fixed covariate model, run by FOCEI
Test6. Use fix() function to introduce a fixed covariate model with a extreme power 100, run by FOCEI