Error: invalid class “omegalist” object: Invalid matrix: determinant is less than 0.

[toothpicky]

Hello,

It does not like my OMEGABLOCK and I do not know why.

$CMT CENTRAL PERIPH CELLS FLU CY
EFFECT

$PARAM

BASET=1000

TVCLNSAT=0.0223076 TVV1=5.40434 TVCLRA=1.43493 TVV2=5.32937 TVCLSAT=2.03386 TVC50=8320.18

TVKOUT=0.0236945 TVEMAX=154.103 TVIC50=895.295 TVGAIN=0.00842024 TVPOWER=2.26102

TVFLMAX=2.49234 TVCYMAX=1.34301

KE0=11.8418

FC50=2.99573 FPOWER=1.32373 KFLU=3.0204 CC50=0.552306 CPOWER=0.905384 KCY=3.15686

LAG=0.210499

TRUNC=0

$MAIN

double CLNSAT=TVCLNSAT exp(ETACLNSAT) ; double V1=TVV1 exp(ETAV1) ; double CLRA=TVCLRA exp(ETACLRA); double V2=TVV2 exp(ETAV2); double CLSAT=TVCLSAT exp(ETACLSAT) ; double C50=TVC50 exp(ETAC50);

double KOUT=TVKOUT exp(ETAKOUT) ; double BASE=BASET exp(ETABASE) ; double EMAX=TVEMAX exp(ETAEMAX) ; double IC50=TVIC50 exp(ETAIC50) ; double GAIN=TVGAIN exp(ETAGAIN) ; double POWER=TVPOWER exp(ETAPOWER) ;

double FLMAX=TVFLMAX exp(ETAFLMAX) ; double CYMAX=TVCYMAX exp(ETACYMAX) ;

ALAG_FLU=LAG; ALAG_CY=LAG; F_CENTRAL=1; F_CELLS=BASE;

$GLOBAL

define LOGT (log10(BASE))

define KNSAT (CLNSAT/ V1)

define KSAT (CLSAT / V1)

define K12 (CLRA / V1)

define K21 (CLRA / V2)

define CP (CENTRAL / V1)

define CE (PERIPH/ V1)

define UNINHIB (1/ (1 + CE / C50))

define KSUM (KNSAT + KSAT * UNINHIB)

define ELIM (CENTRAL * KSUM)

define DIST (CENTRAL K12 - PERIPH K21)

define AFXEFF (1 / (1 + pow(IC50/CE,POWER)))

define AEFFECT (1 + EMAX * AFXEFF)

define FFXEFF (1 / (1 + pow(FC50 / FLU,FPOWER)))

define FEFFECT (1 + FLMAX * FFXEFF)

define CFXEFF (1 / (1 + pow(CC50 / CY,CPOWER)))

define CEFFECT (1 + CYMAX * CFXEFF)

define EFFECT (AEFFECT FEFFECT CEFFECT)

define BASEINP (KOUT * BASE)

define OUT (KOUT CELLS EFFECT)

define IN (BASEINP * pow(BASE /CELLS,GAIN))

if (CELLS<0) CELLS=1; if (CE<0) AFXEFF=0; if (FLU<0.01) FFXEFF=0; if (CY<0.01) CFXEFF=0; if (EFFECT<0) { EFFECT=0; TRUNC=1; }

$ODE

dxdt_CENTRAL = -ELIM - DIST; dxdt_PERIPH = DIST; dxdt_CELLS=IN- OUT; dxdt_FLU=-KFLUFLU; dxdt_CY=-KCYCY; dxdt_EFFECT=KE0CENTRAL - KE0EFFECT ;

$OMEGA @block @labels ETACLNSAT ETAV1 ETACLRA ETAV2 ETACLSAT ETAC50 1.39024 -0.893117 0.573763 0.257089 -0.162566 0.878871 -0.100775 0.0658327 0.346015 0.195481 -0.834404 0.536852 0.0973992 0.196067 1.12526 0.879892 -0.565874 0.0845373 -0.183245 -1.39345 2.78636

$OMEGA @block @labels ETAKOUT ETABASE ETAEMAX ETAIC50 ETAGAIN 1.92093 0.0924729 0.0279971 -0.564586 -0.0440894 1.49413 -0.00249092 -0.0622271 1.17805 1.84527 0.0335413 0.00771135 -0.00928018 -0.0225095 0.305433

$OMEGA @labels ETAFLMAX ETACYMAX 1.02803 0.42464

$SIGMA @labels PROP ADD 0.234989 0.0115461

[FelicienLL]

Hi @toothpicky, Your first matrix is not positive definite.

library(mrgsolve)
mat <- bmat(
  1.39024, 
  -0.893117, 0.573763,
  0.257089, -0.162566, 0.878871,
  -0.100775, 0.0658327, 0.346015, 0.195481,
  -0.834404, 0.536852, 0.0973992, 0.196067, 1.12526,
  0.879892, -0.565874, 0.0845373, -0.183245, -1.39345, 2.78636
)
det(mat)
#> [1] -1.387182e-08

What version of mrgsolve are you using? Up to version 1.0.4, mrgsolve did not handle these small yet negative values, as you can see below.

packageVersion("mrgsolve")
#> [1] '1.0.3'
mrgsolve::mvgauss(mat)
#>       [,1] [,2] [,3] [,4] [,5] [,6]
#>  [1,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [2,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [3,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [4,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [5,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [6,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [7,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [8,]  NaN  NaN  NaN  NaN  NaN  NaN
#>  [9,]  NaN  NaN  NaN  NaN  NaN  NaN
#> [10,]  NaN  NaN  NaN  NaN  NaN  NaN

Maybe it can be fixed by updating the package for a more recent version.

packageVersion("mrgsolve")
#> [1] '1.0.6'
mrgsolve::mvgauss(mat)
#>             [,1]        [,2]       [,3]        [,4]       [,5]       [,6]
#>  [1,] -2.2902972  1.46830525 -1.4065275 -0.42721212  0.9103399 -0.3909998
#>  [2,]  0.1406025 -0.09254095 -0.3984435 -0.03255888 -0.3509535  1.3047007
#>  [3,]  0.8458820 -0.53958196  1.2986216  0.29757575  0.9712909 -0.7168624
#>  [4,] -1.3095470  0.84041994 -0.5145891  0.12888058  0.2977376 -1.1183475
#>  [5,]  1.1884248 -0.77011826 -2.1927070 -1.37993793 -0.8670847 -0.2946304
#>  [6,]  0.2717729 -0.17446648  0.1277216 -0.01120173  0.7943491 -0.5560711
#>  [7,]  1.9901197 -1.28024324 -0.1600081 -0.65553009 -1.9854490  3.7516758
#>  [8,] -1.5046170  0.96709019 -0.1743852  0.27545693  2.8031851 -4.1622630
#>  [9,]  0.7460116 -0.47557465  1.3029838  0.02263070 -0.2326658  2.7258805
#> [10,] -1.0408432  0.67443923  1.8725868  1.22646498  1.1784742 -0.6751965

Félicien

[toothpicky]

I just downloaded mgrsolve yesterday so I am using the latest version. Any workarounds? It runs fine in NONMEM.

packageVersion("mrgsolve") ‘1.0.6’

[kylebaron]

Hi -

At the moment, mrgsolve won't take a matrix that has a determinant less than zero. One thing you could do is coerce the matrix so that it is positive definite. If the matrix is coming from NONMEM, you can use $NMEXT block to import the matrix right from the NONMEM results; only very rarely will that not work and it's usually with a NONMEM run that ends with problems.

Kyle

library(mrgsolve)
#> 
#> Attaching package: 'mrgsolve'
#> The following object is masked from 'package:stats':
#> 
#>     filter
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union

mat <- bmat(
1.39024,
-0.893117, 0.573763,
0.257089 ,-0.162566 ,0.878871,
-0.100775, 0.0658327, 0.346015, 0.195481,
-0.834404, 0.536852, 0.0973992, 0.196067, 1.12526,
0.879892, -0.565874, 0.0845373 ,-0.183245, -1.39345, 2.78636
)
mat <- Matrix::nearPD(mat)$mat
cat(mat[lower.tri(mat, diag=TRUE),]) 
#> 1.39024 -0.8931166 0.257089 -0.100775 -0.834404 0.879892 0.5737635 -0.162566 0.0658327 0.536852 -0.565874 0.878871 0.346015 0.0973992 0.0845373 0.195481 0.196067 -0.183245 1.12526 -1.39345 2.78636
det(as.matrix(mat))
#> [1] 9.111982e-10
eigen(as.matrix(mat))$values
#> [1] 4.419471e+00 1.366168e+00 9.700533e-01 1.740597e-01 2.022435e-02
#> [6] 4.419470e-08

^{Created on 2022-09-08 with reprex v2.0.2}

[toothpicky]

I wrote the simulation file in NONMEM and it runs fine. I assume $NMEXT reads from the NONMEM model outputs. It's not that easy since there's two models, a PK and PD models where I got my THETA, OMEGA, and SIGMAs. I'm C&P the OMEGA straight from my NONMEM outputs so I am not sure why it's not working. I even tried OMEGAs from a different model run and got the same error message.

Coercing it into positive works, I will check the simulations with the NONMEM runs to see how the results compare.

Forgive me, I'm brand new to mgrsolve but now it's showing a different error message:

mod<-mread('mrgsolve/TC') Building TC ... error.

---:: stderr ::--------------------------------------------- Warning message: In system(cmd) : 'make' not found

Error: the model build step failed.

Do you see any glaring errors in this model?

$CMT CENTRAL PERIPH CELLS FLU CY
EFFECT

$PARAM

BASET=1000

TVCLNSAT=0.0223076 TVV1=5.40434 TVCLRA=1.43493 TVV2=5.32937 TVCLSAT=2.03386 TVC50=8320.18

TVKOUT=0.0236945 TVEMAX=154.103 TVIC50=895.295 TVGAIN=0.00842024 POWER=2.26102

TVFLMAX=2.49234 TVCYMAX=1.34301

KE0=11.8418

FC50=2.99573 FPOWER=1.32373 KFLU=3.0204 CC50=0.552306 CPOWER=0.905384 KCY=3.15686

LAG=0.210499

TRUNC=0

$MAIN

double CLNSAT=TVCLNSAT exp(ETACLNSAT) ; double V1=TVV1 exp(ETAV1) ; double CLRA=TVCLRA exp(ETACLRA); double V2=TVV2 exp(ETAV2); double CLSAT=TVCLSAT exp(ETACLSAT) ; double C50=TVC50 exp(ETAC50);

double KOUT=TVKOUT exp(ETAKOUT) ; double BASE=BASET exp(ETABASE) ; double EMAX=TVEMAX exp(ETAEMAX) ; double IC50=TVIC50 exp(ETAIC50) ; double GAIN=TVGAIN * exp(ETAGAIN) ;

double FLMAX=TVFLMAX exp(ETAFLMAX) ; double CYMAX=TVCYMAX exp(ETACYMAX) ;

ALAG_FLU=LAG; ALAG_CY=LAG; F_CENTRAL=1; F_CELLS=BASE;

$GLOBAL

define LOGT (log10(BASE))

define KNSAT (CLNSAT/ V1)

define KSAT (CLSAT / V1)

define K12 (CLRA / V1)

define K21 (CLRA / V2)

define CP (CENTRAL / V1)

define CE (EFFECT/ V1)

define UNINHIB (1/ (1 + CE / C50))

define KSUM (KNSAT + KSAT * UNINHIB)

define ELIM (CENTRAL * KSUM)

define DIST (CENTRAL K12 - PERIPH K21)

define AFXEFF (1 / (1 + pow(IC50/CE,POWER)))

define AEFFECT (1 + EMAX * AFXEFF)

define FFXEFF (1 / (1 + pow(FC50 / FLU,FPOWER)))

define FEFFECT (1 + FLMAX * FFXEFF)

define CFXEFF (1 / (1 + pow(CC50 / CY,CPOWER)))

define CEFFECT (1 + CYMAX * CFXEFF)

define EFFECT (AEFFECT FEFFECT CEFFECT)

define BASEINP (KOUT * BASE)

define OUT (KOUT CELLS EFFECT)

define IN (BASEINP * pow(BASE /CELLS,GAIN))

if (CELLS<0) CELLS=1; if (CE<0) AFXEFF=0; if (FLU<0.01) FFXEFF=0; if (CY<0.01) CFXEFF=0; if (EFFECT<0) { EFFECT=0; TRUNC=1; }

$ODE

dxdt_CENTRAL = -ELIM - DIST; dxdt_PERIPH = DIST; dxdt_CELLS=IN- OUT; dxdt_FLU=-KFLU FLU; dxdt_CY=-KCY CY; dxdt_EFFECT=KE0 CENTRAL - KE0 EFFECT ;

$OMEGA @block @labels ETACLNSAT ETAV1 ETACLRA ETAV2 ETACLSAT ETAC50 1.39024 -0.8931166 0.257089 -0.100775 -0.834404 0.879892 0.5737635 -0.162566 0.0658327 0.536852 -0.565874 0.878871 0.346015 0.0973992 0.0845373 0.195481 0.196067 -0.183245 1.12526 -1.39345 2.78636

$OMEGA @block @labels ETAKOUT ETABASE ETAEMAX ETAIC50 ETAGAIN 1.92093 0.0924729 0.0279971 -0.564586 -0.0440894 1.49413 -0.00249092 -0.0622271 1.17805 1.84527 0.0335413 0.00771135 -0.00928018 -0.0225095 0.305433

$OMEGA @labels ETAFLMAX ETACYMAX 1.02803 0.42464

$SIGMA @labels PROP ADD 0.234989 0.0115461

[kylebaron]

@toothpicky -

I am not sure why it's not working

It's not working because the determinant of the first OMEGA matrix is negative; it also not positive definite; mrgsolve checks this and asks that you put in positive definite matrix for safety. 99.999% of the time if you import directly from NONMEM it will give you a well formed matrix.

Where are you copying the values from? You probably just need a couple more digits on each number to get it to work.

You can import from two different models ... PK and PD ... into the same mrgsolve control stream. I'd encourage you to go that route. But let's look at the .ext file first. Can you share that from the first model?

It appears that you are getting the new error message because you don't have Rtools installed? Assuming you are on Windows computer. Did you do this?

Depending on what version of R you have, install from here: https://cran.r-project.org/bin/windows/Rtools/

Let me know if you have any issues.

Kyle

[toothpicky]

Thanks Kyle. I was using R tools 4.2 but my R was version 4.1 so I updated R and the error message went away.

Now there is something wrong with my if statements?

Building TC ... error.

---:: stderr ::--------------------------------------------- In file included from TC-mread-source.cpp:3: TC-mread-header.h:104: warning: "EFFECT" redefined 104	#define EFFECT A[5]

TC-mread-header.h:58: note: this is the location of the previous definition 58	#define EFFECT (AEFFECT FEFFECT CEFFECT)

28:20: error: assignment of read-only location '(& A)->std::vector::operator' 28 | if (CELLS<0) {CELLS=1;} In file included from TC-mread-source.cpp:3: TC-mread-header.h:52:19: error: lvalue required as left operand of assignment 52 | #define AFXEFF (1 / (1 + pow(IC50/CE,POWER))) | ~^~~~~~~~~ 29:11: note: in expansion of macro 'AFXEFF' 29 | if (CE<0) AFXEFF=0; | ^~ TC-mread-header.h:54:19: error: lvalue required as left operand of assignment 54 | #define FFXEFF (1 / (1 + pow(FC50 / FLU,FPOWER))) | ~^~~~~~~~~ 30:15: note: in expansion of macro 'FFXEFF' 30 | if (FLU<0.01) FFXEFF=0; | ^~ TC-mread-header.h:56:19: error: lvalue required as left operand of assignment 56 | #define CFXEFF (1 / (1 + pow(CC50 / CY,CPOWER))) | ~^~~~~~~~ 31:15: note: in expansion of macro 'CFXEFF' 31 | if (CY<0.01) CFXEFF=0; | ^~ 33:8: error: assignment of read-only location '(& A)->std::vector::operator' 33 | EFFECT=0; 34:7: error: assignment of read-only location '(& THETA)->std::vector::operator' 34 | TRUNC=1; In file included from TC-mread-header.h:14, from TC-mread-source.cpp:3: At global scope: C:/Users/anh.nguyen/AppData/Local/R/win-library/4.2/mrgsolve/base/modelheader.h:68:19: error: expected declaration before '}' token 68 | #define DONE }} | ^ C:/Users/anh.nguyen/AppData/Local/R/win-library/4.2/mrgsolve/base/modelheader.h:65:22: note: in expansion of macro 'DONE' 65 | #define END_main DONE | ^~~~ 37:1: note: in expansion of macro 'END_main' 37 | __END_main__ | ^~~~ make: *** [C:/PROGRA~1/R/R-42~1.1/etc/x64/Makeconf:260: TC-mread-source.o] Error 1

Error: the model build step failed.

[kylebaron]

Some issues:

FIrst

if (CELLS<0) {CELLS=1;}

You can't re-assign the value of a compartment. If you want you can say

double CELLS2 = CELLS;
if(CELLS < 0) CELLS2 = 1;

and then use CELLS2 in your model in place of CELLS.

Second

Pre-processor definitions are one-way; you're just aliasing some name to a bit of code. You shouldn't assign anything to it.

#define AFXEFF (1 / (1 + pow(IC50/CE,POWER)))

...

 if (CE<0) AFXEFF=0;

Third

You can't really write this under global (beyond re-assigning compartments and pp directives); I'd ailas these with another variable (like I showed above). Depending on what you're trying to do with it, it could be in $TABLE or $MAIN.

$GLOBAL 
...
if (CELLS<0) CELLS=1;
if (CE<0) AFXEFF=0;
if (FLU<0.01) FFXEFF=0;
if (CY<0.01) CFXEFF=0;
if (EFFECT<0) {
EFFECT=0;
TRUNC=1;
}

I think this can all get fixed; it's not that far. I can look at it with you to better understand the intent and get the code lined up. If you have a NONMEM model to work from that would be easiest. If you want to email and do it off line, we could do that too.

Kyle

[toothpicky]

Thanks Kyle! I have sent an email to you at 'kylebtwin@imap.cc'

[toothpicky]

Hello, I am still having an issue with negative determinants. I turned it into a positive determinant and it's still claiming it's less than 0.

mat <- bmat( 2.50E+00, 1.09E+00 , 5.03E-01, -6.03E-01, -2.06E-01, 2.62E-01, -4.23E-01, -1.58E-01, 1.56E-01, 1.03E-01, -8.57E-01, -3.38E-01, 2.84E-01, 1.21E-01, 1.06E+00, 1.00E+00, 4.16E-01, -2.91E-01, -1.63E-02, -2.40E+00, 6.19E+00 ) det(mat)

-4.623566e-10

mat <- Matrix::nearPD(mat)$mat cat(mat[lower.tri(mat, diag=TRUE),])

2.50002 1.089964 -0.602985 -0.4229835 -0.8570107 0.9999958 0.5030634 -0.2060268 -0.1580294 -0.3379809 0.4160075 0.2620113 0.1560124 0.283992 -0.2910032 0.1030136 0.1209912 -0.01630349 1.060006 -2.399998 6.190001

det(as.matrix(mat))

3.051959e-13

[toothpicky]

My matrix can be read by this mvgauss function but there is still an error message for negative determinants, even when the determinant is positive.

 mrgsolve::mvgauss(mat)