Closed billdenney closed 2 years ago
billdenney
commented
2 years ago My preference for the behavior would be that the user would be required to set an option to either initialize at 0 or at the first measurement. And, that setting should be an explicit choice rather than a default.
mattfidler
commented
2 years ago There is already an option, though perhaps documentation needs to be a bit more explicit. It is very explicit when the first observation is negative,
https://github.com/nlmixr2/rxode2/blob/6d440074b4eefe5adce79887739f59bf2d4932e4/src/etTran.cpp#L2045
So use:
rxSetIni0(FALSE)To change this behavior
mattfidler
commented
2 years ago Evetually, I want to change to something like:
a(10)=3
b(4)=6But that has been a feature request that I have been thinking about for awhile (but never created an issue for)...
mattfidler
commented
2 years ago Because the option is already present, there is likely a data translation error or there is an error for the initial problem somehow
mattfidler
commented
2 years ago Well it isn't a data translation error. It is doing the right thing (using the underlying classic rxode2 code and the prepfit data.
Note that I still can't run reprex::reprex() on wayland:
> library(rxode2)
> # from (oncology_sdm_lobo_2002())$simulationModel
+ rx <- rxode2({
+ param(lkng, ltau, lec50, kmax, propErr, addErr, cp, tumor0)
+ kng = exp(lkng)
+ tau = exp(ltau)
+ taulast = tau
+ ec50 = exp(lec50)
+ edrug = kmax * cp/(ec50 + cp)
+ tumor(0) = tumor0
+ d/dt(transit1) = (edrug - transit1)/tau
+ d/dt(transit2) = (transit1 - transit2)/tau
+ d/dt(transit3) = (transit2 - transit3)/tau
+ d/dt(transitlast) = transit3/tau - transitlast/taulast
+ d/dt(tumor) = kng * tumor - transitlast * tumor
+ rx_yj_ ~ 2
+ rx_lambda_ ~ 1
+ rx_low_ ~ 0
+ rx_hi_ ~ 1
+ rx_pred_f_ ~ tumor
+ rx_pred_ ~ rx_pred_f_
+ rx_r_ ~ (addErr)^2 + (rx_pred_f_)^2 * (propErr)^2
+ ipredSim = rxTBSi(rx_pred_, rx_lambda_, rx_yj_, rx_low_,
+ rx_hi_)
+ sim = rxTBSi(rx_pred_ + sqrt(rx_r_) * err.tumor, rx_lambda_,
+ rx_yj_, rx_low_, rx_hi_)
+ dvid(5)
+ })
+ > # from qs::base91_encode(qs::qserialize(d_sim))
+ prepfit <- qs::qdeserialize(qs::base91_decode("un]\"BAAA@QRtHACAAAAAAAuWeBAABdk1kus^^d8Ah9}?=Z:alBMc4Iv(F\":C?hVBAMZRxwFfBB7IB.y6FTL)yFQA@v;KlgSH3Vn~rL/,{CP/ez~`.3>$Wj$rcy==//#}Pu?\"V(RgKtf1J3qQg/yI7*1]/WV=iegVmPs3?a\":kEMu~/*zmX#;E4`i@It`]Ouu[N]T8G3!4A.^j0<Zd;mY7D`P6Z|KTj_p?r7u[IPuyFRoUcL\"6u|(G_on6g9c{ZLJ[_gE^&47rbL(#6W{EA%hQUW!][2;k<}\"(4SB{!~>M%xmxwMU%ET6#~xvX9rH!;S53gbLTWY]Fcri\"]7\"|Z^W{xobiiTc~DLN_;.Itj(INGKCupDYxEA^!GzfHO=aDW&(I)z}0*mZD\"^b.O!QdY2rVRD;~Z*HB(]G_Dpj]*0A`]+7VoGDF@,vk>jx}tFI>MVOnZojuABN9Bt\"O~V[n6U[kn|W74&xR7CL(Skn:CA)NP`||hQ%w/i+&c8$#KxsFdb4,qI\"Fl&lLg,?$eh&s{`QxtwPWi$GX<[*<0{to@[:NAy}a=O`wedEA*Abqhz2bL2sfII3ZRJR#5q~:FeBW%/F<]`(?Q:c(qc,DZ_d.&|J(NW~Q4kz;Us(7e+Z0YGMdvf.%XRgD]FA2D10sl^KxuPXvXSm+p}ndVY!3`o}Iq+M;i~mLmr1In0~ymm]K2x9g9Ij.UkBOTriq+93<po9tNj@%%W#FzA+/MMb]k3YmXS*RdjE{?pjnr%q6}&81.2&#ni.Au{pL>#eKT9uaMmvF7L^aDwL?sj>s|}[XMQdEx(yS|vIDfwqH:YXc(EfrzuplGn3`|X=ObNnD%;3(ST3tWr^D+vDG=cjKk!^:5ZfpXK5/dqF@dW6+*lb~@\"*H_t@3rRG9w|kG1!SRghm{sUOcDQ?.gYd?c;:IC~RF6lEfd4X;I3+4Y&8.x0P[h8]Ffo)$Y=7)|08Vid5@2btxd4I[0>E7~+CCX~|{Ve0qf^aQ1M4vsQN<B6]~R{LE*vH*8[Q|b[`QO%/.,mr9d<+O/PqesyX).kR)a[.s^OsSJLRZgrwt?NA#EQf)p,+wjGm!8aRz9^SKrW5`PTsjH_EVK.eZ)T]k0Md!mRLU%+}K!{BB>[uCto@(_b&&NA0eV*R}~Tp<Ey:7|>#Kz+AL~%_Wu@J&>C__{[#o3Lmxh$%R;264SAA)O;,:*:p.s:z+=pot[NkrFE@[}VoBF=P<AI)v_tBGUZy4^}6MUL|glMtnz]K/D|qUE``3x])hpz{k$OXX}C4B%m,AwVR)[N\"d$i}eLLJ;bFzpavPm|k^Ef;$C|S+GBP;CdOQSE<>v:;Cj#YR]M@a{68?;s&4Y4Zd[D`znL|kp1n)}0+Q`(/LcQM42JIA}ZFH+9>#dQLs>JX?#GT<W+#YOo[zpESq=N.@3/bpfI)5wKpkIRY>_M&#wui|N/HWs/WyJ/CSn*BaX!:]**?VL`zPu?xA|.+kn)M[xCjLHQ@8A6HX^T4hXT7F~JQw])fR/^Y>YBD([_QsK;[iYPcbW2Z\"BA9q~@t^rOgPg95OnA01xfA00E8k)eQtyH#MzpQ$4|et<eI&o|Uhx|t19ZFY$DAw5Z<,6p7MQM*N07@b$w):77pmoP{iZ4KvtvMyKYw|?{al}MZNtX<:RTKFB?PAdg=76:~4(HH*s62mDm2K:Qw^(9,*ICI.`?HiTsl8KIYV0C!0%0S:JM6(#2A)w(y(<9\"?n`k{0V!yrmQChNo8WWL=/*eAA_+&bKol/+*eiJi*e2dIW!hi.)+RHLR9O}Z69<<eE;aD%#B0ql|30pbBs3Z^0(q$nWfGs|LD$erJR2r@>6W(eGPRmx0B?bGrM@DA31ha^UZWMe52Egm6nn3=lcs3[gq!vBvrKiTSf],3$wjk7^`BPE#{kPh.CcbM*h:9+XPRfGx1P^)X!utw{#;**|LZtNdAwx\"IC?m[]Wil9!m3Tx}{Xx@9;`=+ODLPp)w]zc\"D.rl7/d:ismtgiesu]/CF[9v7ICZ:U?wky`@wQn])q@|Dd#IfQbwAQF`7]{f;hJCK=Vv)N6M4x4>G.u%o7I5fCay=\"^,?M!600bO\"MBTL9~4}eDkY:YkGc6G9oMzo3<9[Ye`GFvPyXiw5$Z]*PByzRza9ik5LJ;}@p;JH}+.,1om5gxAkHdn%`#kbKMJuiDOa#*M!h.2dUje501>,MIW$207SwS$M~FCMaG:`&l?j<c+ZMyFt+u:CpBWOZL)$@WywmY|LCGPdVv=>K]e\"!>HG]Tw0Xbm7SOIjt.B%cH5Ewy_YgfI6(#$I]1Q4Jte1:^W%XdqR4C#OV#._I`G$k|#>NfxcZp!Y]YnX~yZ0nx),%jfE\"u7yS,:y&oJV`;4nMg0g[a]Tidjt|y[*~_8H;*X61T/1kv9x$f/3N?~>xL$mO^0Nk6W|U[wD"))
> head(prepfit) # time 0 is missing here
id time kng tau taulast ec50 edrug DV sim transit1
1 1 1 0.02 34.1 34.1 50 7.303850e-05 102.0202 65.673434 5.666220e-06
2 1 2 0.02 34.1 34.1 50 2.699399e-05 104.0811 104.740533 7.603913e-06
3 1 3 0.02 34.1 34.1 50 9.947509e-06 106.1837 37.562447 8.164464e-06
4 1 4 0.02 34.1 34.1 50 3.661786e-06 108.3287 3.156474 8.222402e-06
5 1 5 0.02 34.1 34.1 50 1.347408e-06 110.5171 47.265125 8.080313e-06
6 1 6 0.02 34.1 34.1 50 4.957259e-07 112.7497 117.696059 7.880630e-06
transit2 transit3 transitlast tumor cp tumor0
1 8.225240e-08 1.011216e-09 1.243437e-11 102.0202 0.367879441 100
2 2.727191e-07 6.000752e-09 1.006670e-10 104.0811 0.135335283 100
3 4.928050e-07 1.688075e-08 4.147003e-10 106.1837 0.049787068 100
4 7.159394e-07 3.388743e-08 1.123073e-09 108.3287 0.018315639 100
5 9.308434e-07 5.674403e-08 2.388233e-09 110.5171 0.006737947 100
6 1.134605e-06 8.498855e-08 4.356776e-09 112.7497 0.002478752 100
> head(etTrans(prepfit, rx)) # time 0 is added back here with internal evid=9, the evid to initialize to zero.
ID TIME EVID AMT II DV cp
1 1 0 9 NA 0 NA NA
2 1 1 0 NA 0 102.0202 0.367879441
3 1 2 0 NA 0 104.0811 0.135335283
4 1 3 0 NA 0 106.1837 0.049787068
5 1 4 0 NA 0 108.3287 0.018315639
6 1 5 0 NA 0 110.5171 0.006737947
Covariates (non time-varying):
ID tumor0
1 1 100
2 2 100
Compartment translation:
[1] Compartment Number
<0 rows> (or 0-length row.names)
> Note that turning off the ini0, the model runs fine:
> library(tidyverse)
+ library(nlmixr2)
── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.2 ──
✔ ggplot2 3.3.6 ✔ purrr 0.3.5
✔ tibble 3.1.8 ✔ dplyr 1.0.10
✔ tidyr 1.2.1 ✔ stringr 1.4.1
✔ readr 2.1.3 ✔ forcats 0.5.2
── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
> Loading required package: nlmixr2data
> . + . + > d_sim <-
+ data.frame(
+ tumor0 = 100,
+ time = 0:24,
+ cmt = "tumor",
+ cp=pmxTools::calc_sd_1cmt_linear_bolus(t = 0:24, dose = 1, CL = 1, V = 1)
+ ) %>%
+ crossing(
+ ID = 1:2
+ )
+ > simulated <- nlmixr2(oncology_sdm_lobo_2002, data = d_sim, est = "rxSolve")
+ prepfit <-
+ as.data.frame(simulated) %>%
+ filter(!(id == 1 & time == 0)) %>%
+ rename(
+ DV=ipredSim
+ )
+ > rxode2::rxSetIni0(FALSE)
[1] FALSE
> simfit <- nlmixr2(oncology_sdm_lobo_2002, data= prepfit, est = "focei", control = foceiControl(print = 0))
rxode2 2.0.9 using 4 threads (see ?getRxThreads)
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
✔ done
→ compress origData in nlmixr2 object, save 6624
→ compress parHist in nlmixr2 object, save 5296
> simfit
── nlmixr² Population Only (outer: nlminb) ──
OBJF AIC BIC Log-likelihood Condition Number
Pop 62.79619 164.8522 176.2031 -76.42608 80745.15
── Time (sec simfit$time): ──
setup optimize covariance table compress other
elapsed 0.016787 0.006418 0.00642 0.021 0.011 0.197375
── (simfit$parFixed or simfit$parFixedDf): ──
Parameter Est.
lkng Cell net growth rate (growth minus death) (1/hr) -3.88
ltau Mean transit time of each transit compartment (hr) 3.56
lec50 Drug concentration reducing the cell growth by 50% (ug/mL) 3.92
kmax Maximum drug-related reduction in cell growth (1/hr) -19
propErr Proportional residual error (fraction) 0.0115
addErr Additive residual error (tumor volume units) 0.0973
SE %RSE Back-transformed(95%CI) BSV(SD) Shrink(SD)%
lkng 0.0137 0.353 0.0206 (0.0201, 0.0212)
ltau 0.27 7.58 35.2 (20.8, 59.8)
lec50 0.257 6.55 50.2 (30.4, 83.1)
kmax 3.87 20.4 -19 (-26.5, -11.4)
propErr 0.0115
addErr 0.0973
Covariance Type (simfit$covMethod): |r|
Information about run found (simfit$runInfo):
• IDs without zero-time start at the first observed time: 1
• gradient problems with initial estimate and covariance; see $scaleInfo
• R matrix non-positive definite but corrected by R = sqrtm(R%*%R)
Censoring (simfit$censInformation): No censoring
Minimization message (simfit$message):
false convergence (8)
In an ODE system, false convergence may mean "useless" evaluations were performed.
See https://tinyurl.com/yyrrwkce
It could also mean the convergence is poor, check results before accepting fit
You may also try a good derivative free optimization:
nlmixr2(...,control=list(outerOpt="bobyqa"))
── Fit Data (object simfit is a modified tibble): ──
# A tibble: 49 × 16
ID TIME DV IPRED IRES IWRES transit1 transit2 transit3 transi…¹ tumor
<fct> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 1 1 102. 100 2.02 1.76 0 0 0 0 100
2 1 2 104. 102. 2.00 1.70 -0.00386 -5.45e-5 -5.27e-7 -3.98e-9 102.
3 1 3 106. 104. 1.98 1.65 -0.00517 -1.80e-4 -3.72e-6 -5.56e-8 104.
# … with 46 more rows, 5 more variables: kng <dbl>, tau <dbl>, taulast <dbl>,
# ec50 <dbl>, edrug <dbl>, and abbreviated variable name ¹transitlast
# ℹ Use `print(n = ...)` to see more rows, and `colnames()` to see all variable names
> For me when I change the option I get:
> simfit <- nlmixr2(oncology_sdm_lobo_2002, data= prepfit, est = "focei", control = foceiControl(print = 0))
calculating covariance matrix
[====|====|====|====|====|====|====|====|====|====] 100%; 0:00:00 done
→ Calculating residuals/tables
DLSODA- At start of problem, too much accuracy
requested for precision of machine.. See TOLSF (=R1)
IDID=-3, illegal input detected (see printed message).
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 1000000 are needed
IDID=-2, larger nmax is needed
✔ done
→ compress origData in nlmixr2 object, save 6624
→ compress parHist in nlmixr2 object, save 5008Note that
DLSODA- At start of problem, too much accuracy
requested for precision of machine.. See TOLSF (=R1)
IDID=-3, illegal input detected (see printed message).
exit of dop853 at x = 0.0000000000000000e+00, more than nmax = 1000000 are needed
Says that this is an ODE solving issue in general at the table step, so it is related to nlmixr2/nlmixr2est#256 instead.
mattfidler
commented
2 years ago Or rather that may be the cause.
mattfidler
commented
2 years ago As suggested by the initial issue, the model had the issue when removing the subject's evid=9. My guess is that the estimation no longer treated the first data point (because it no longer is part of the data) and had different estimates, and then had problems solving the resulting ode
billdenney
commented
2 years ago FYI, all of my recent issues are related to this problem. I think that fixing the data handling around time = 0 and time != 0 will resolve my issues for all of these. And for now, I can work around it with an EVID = 2 row at the beginning.
mattfidler
commented
2 years ago Ok. I tried to reproduce it again and I couldn't. I am reinstalling the entire stack to see if I can reproduce again.
mattfidler
commented
2 years ago With evid=9
Browse[2]> etTrans(fit$dataSav, model)
ID TIME EVID AMT II DV cp tumor0
1 1 0 9 NA 0 NA NA NA
2 1 1 0 NA 0 102.0202 3.678794e-01 100
3 1 2 0 NA 0 104.0811 1.353353e-01 100
4 1 3 0 NA 0 106.1837 4.978707e-02 100
5 1 4 0 NA 0 108.3287 1.831564e-02 100
6 1 5 0 NA 0 110.5171 6.737947e-03 100with evid=0
Browse[2]> etTrans(fit$dataSav, model)
ID TIME EVID AMT II DV cp
1 1 0 0 NA 0 100.0000 1.000000e+00
2 1 1 0 NA 0 102.0202 3.678794e-01
3 1 2 0 NA 0 104.0811 1.353353e-01
4 1 3 0 NA 0 106.1837 4.978707e-02
5 1 4 0 NA 0 108.3287 1.831564e-02
6 1 5 0 NA 0 110.5171 6.737947e-03
7 1 6 0 NA 0 112.7497 2.478752e-03
8 1 7 0 NA 0 115.0274 9.118820e-04
9 1 8 0 NA 0 117.3512 3.354626e-04
...
Covariates (non time-varying):
ID tumor0
1 1 100
2 2 100with evid=2
etTrans(fit$dataSav, model)
ID TIME EVID AMT II DV cp
1 1 0 101 0 0 NA 1.000000e+00
2 1 0 2 NA 0 NA 1.000000e+00
3 1 1 0 NA 0 102.0202 3.678794e-01
4 1 2 0 NA 0 104.0811 1.353353e-01
5 1 3 0 NA 0 106.1837 4.978707e-02
6 1 4 0 NA 0 108.3287 1.831564e-02
7 1 5 0 NA 0 110.5171 6.737947e-03
8 1 6 0 NA 0 112.7497 2.478752e-03
9 1 7
...
Covariates (non time-varying):
ID tumor0
1 1 100
2 2 100
Compartment translation:
[1] Compartment Number
<0 rows> (or 0-length row.names)Which means imputed variables need to be considered before observed.
mattfidler
commented
2 years ago They are already imputed. The only issue is time invariant covariates are not considered correctly with NA values
In the example below, the first subject in
simfithas no simulated results because the time 0 result is removed. (I think that this is not an issue when the first time is not zero for both subjects.)Created on 2022-10-17 with reprex v2.0.2