There is error in simulation with covariance matrix, I don't know what's the problem is. #972

WXX155367 · 2025-12-17T12:36:10Z

WXX155367
Dec 17, 2025

Hi guys. I am doing simulation with covariance matrix. I followed tutorial (11.3.2) and wrote "omega" in "theta" and also wrote covariance of omega in "thetaMat" because I want to consider the variance during simulation. But it showed the error as below:
"warning: chol(): given matrix is not symmetric
error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits, :
error calculating 'chol(thetaMat)'
I would like to know if there is any mistake in my code or not. I put my code as below:

#build PK model
library(rxode2)
library(dplyr)
#library(posologyr)
rx1 <- rxode2({
    TVCL <- tvcl*(EGFR/52.0)^tvcl_egfr*(BMI/24.0)^tvcl_bmi
    CL <- TVCL*exp(eta.cl)
    TVV1 <- tvv1*(BMI/24.0)^tvv1_bmi*(TOPR/66.0)^tvv1_topr*(SODI/138.25)^tvv1_sodi*(1+tvv1_fb*(FB-2.7965))*tvv1_flg
    V1 <- TVV1*exp(eta.v1)
    Q <- tvq
    V2 <- tvv2
    A1 ~ centr;
    A2 ~ peri;
    d/dt(centr) <- -A1*(CL/V1 + Q/V1) + A2*Q/V2;
    d/dt(peri) <- A1*Q/V1 - A2*Q/V2;
    ipred <- centr/V1;
    Dv <- ipred + ipred*prop.sd*eps1
})

theta <- c(prop.sd=0.189,
           tvcl=1.81,
           tvv1=16,
           tvq=0.769,
           tvv2=12.6,
           tvcl_egfr=0.856,
           tvcl_bmi=0.858,
           tvv1_bmi=0.421,
           tvv1_fb=0.0000451,
           tvv1_flg=0.903,
           tvv1_topr=-0.529,
           tvv1_sodi=-2.13,
           eta.cl=0.14,
           eta.v1=0.04)

sigma <- lotri(eps1 ~ 1)   # variance=1

#simulated with covariance matrix
thetaMat <- lotri(
  prop.sd + tvcl + tvv1 + tvq + tvv2 + tvcl_egfr + tvcl_bmi + tvv1_bmi + tvv1_fb + tvv1_flg + tvv1_topr + tvv1_sodi + eta.cl + eta.v1~
    c(0.000486,
      -0.000534, 0.0139,
      -0.00258, 0.0225, 0.269,
      0.000845, -0.00527, -0.00809, 0.00751,
      0.0263, -0.418, -0.479, 0.32, 0.208,
      0.00039, -0.00545, -0.00607, 0.00611, 0.335, 0.0108,
      -0.000606, 0.000954, 0.0151, 0.00167, 0.104, 0.00836, 0.0309,
      -0.000145, -0.000984, -0.0166, -0.000275, 0.0548, 0.000295, 0.0113, 0.0284,
      -0.0000000274, 0.00000000391, 0.00000127, -0.000000134, 0.000000296, -0.000000204, 0.0000000143, 0.000000448, 0.000000000128,
      0.000105, -0.000714, -0.0137, -0.000241, 0.00985, -0.000194, -0.00118, 0.00253, -0.0000000584, 0.00511,
      0.000228, 0.000948, 0.00705, -0.0013, -0.0379, -0.000816, -0.00151, -0.00753, -0.000000197, -0.00132, 0.024,
      0.000958, -0.00452, -0.0108, 0.00556, 0.00729, 0.0025, -0.00871, 0.00649, -0.00000207, 0.00492, -0.00363, 0.356,
      0.000396, -0.00313, -0.00614, 0.00148, 0.081, -0.000111, -0.00219, 0.000206, 0.000000026, 0.000287, -0.0000323, 0.00203, 0.00184,
      -0.000215, 0.000461, 0.00443, -0.000661, -0.0237, -0.000494, 0.000357, 0.000208, 0.000000051, -0.0000936, 0.000192, -0.00116, -0.000179, 0.000026
    )
)


#create patients data set
library(tmvtnorm)
set.seed(1234)

## 1. define variables
cov_names <- c("BW",                  # Body weight (kg)
               "BMI",                 # Body mass index (kg/m2)
               "FB",                  # 24h fluid balance (L)
               "EGFR",                # eGFR (mL/min/1.73m2)
               "SODI",                # Serum sodium (mmol/L)
               "TOPR")                # Serum total protein (g/L)

## 2. median, max and min
mu     <- c(70.0,   24.7,   2.8,   53.0,   138.4,  66.0)   
lower  <- c(36.0,   13.7,  -1.5,    6.0,   121.9,  26.0)   
upper  <- c(115.0,  39.6,  13.9,  200.0,   155.9,  90.0)  

## 3. covariates of matrix
Sigma <- matrix(
  c(256.9,   69.4,   4.3,  -153.7,   -9.3,   24.8,
     69.4,   25.1,   1.1,   -60.7,   -1.3,   11.1,
      4.3,    1.1,   5.1,   -18.8,    0.9,    0.9,
   -153.7,  -60.7, -18.8,  1375.8,  -31.4,  -98.4,
     -9.3,   -1.3,   0.9,   -31.4,   40.5,   -1.4,
     24.8,   11.1,   0.9,   -98.4,   -1.4,  123.2),
  nrow = 6, byrow = TRUE
)
colnames(Sigma) <- cov_names
rownames(Sigma) <- cov_names

## 4.extract 200 patients
n_virt <- 200

X_trans <- rtmvnorm(
  n      = n_virt,
  mean   = mu,
  sigma  = Sigma,
  lower  = lower, 
  upper  = upper
)
colnames(X_trans) <- cov_names

## 5. create virtual patients data frame
virtual_patients <- as.data.frame(X_trans)
virtual_patients$id <- 1:n_virt
virtual_patients <- virtual_patients[, c("id", cov_names)]

#dose and event design
s <- seq(0, 24, by = 8) 

#create event table
e <- et() %>%
  et(id=virtual_patients$id, amt=15*virtual_patients$BW, duration=0.5, rate=(15*virtual_patients$BW)/0.5, time=0) %>%
  et(s) %>%
  as.data.frame %>%
  merge(virtual_patients, by="id") %>%
  as_tibble

#create simulation
sim <- rxSolve(
  rx1, theta, e,
  thetaMat = thetaMat,
  thetaLower=0,
  omega = c("eta.cl", "eta.v1"),
  omegaXform = "variance",
  sigma = sigma,
  sigmaXform = "identity",
  nStud = 100,
  dfSub=199,
  dfObs=599
)
print(sim)

Answered by mattfidler

Dec 17, 2025

Here is a reproducible solution:

#build PK model
library(rxode2)
#> rxode2 5.0.1.9000 using 11 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
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
#library(posologyr)

rx1 <- rxode2({
  TVCL <- tvcl*(EGFR/52.0)^tvcl_egfr*(BMI/24.0)^tvcl_bmi
  CL <- TVCL*exp(eta.cl)
  TVV1 <- tvv1*(BMI/24.0)^tvv1_bmi*(TOPR/66.0)^tvv1_topr*(SODI/138.25)^tvv1_sodi*(1+tvv1_fb*(FB-2.7965))*tvv1_flg
  V1 <- TVV1*exp(eta.v1)
  Q <- tvq
  V2 <- tvv2
  A1 ~ centr;

View full answer

WXX155367 · 2025-12-17T13:08:54Z

WXX155367
Dec 17, 2025
Author

It seems like I solved that problem. I add code from previous discussion:

thetaMat1 <- as.matrix(Matrix::nearPD(thetaMat)$mat)
thetaMatChol <- chol(thetaMat1)

and it works and the output has sim$params and simulation table.

0 replies

mattfidler · 2025-12-17T21:58:14Z

mattfidler
Dec 17, 2025
Maintainer

Here is a reproducible solution:

#build PK model
library(rxode2)
#> rxode2 5.0.1.9000 using 11 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
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
#library(posologyr)

rx1 <- rxode2({
  TVCL <- tvcl*(EGFR/52.0)^tvcl_egfr*(BMI/24.0)^tvcl_bmi
  CL <- TVCL*exp(eta.cl)
  TVV1 <- tvv1*(BMI/24.0)^tvv1_bmi*(TOPR/66.0)^tvv1_topr*(SODI/138.25)^tvv1_sodi*(1+tvv1_fb*(FB-2.7965))*tvv1_flg
  V1 <- TVV1*exp(eta.v1)
  Q <- tvq
  V2 <- tvv2
  A1 ~ centr;
  A2 ~ peri;
  d/dt(centr) <- -A1*(CL/V1 + Q/V1) + A2*Q/V2;
  d/dt(peri) <- A1*Q/V1 - A2*Q/V2;
  ipred <- centr/V1;
  Dv <- ipred + ipred*prop.sd*eps1
})

theta <- c(prop.sd=0.189,
           tvcl=1.81,
           tvv1=16,
           tvq=0.769,
           tvv2=12.6,
           tvcl_egfr=0.856,
           tvcl_bmi=0.858,
           tvv1_bmi=0.421,
           tvv1_fb=0.0000451,
           tvv1_flg=0.903,
           tvv1_topr=-0.529,
           tvv1_sodi=-2.13,
           eta.cl=0.14,
           eta.v1=0.04)

sigma <- lotri(eps1 ~ 1)   # variance=1

#simulated with covariance matrix
thetaMat <- lotri(
  prop.sd + tvcl + tvv1 + tvq + tvv2 + tvcl_egfr + tvcl_bmi + tvv1_bmi + tvv1_fb + tvv1_flg + tvv1_topr + tvv1_sodi + eta.cl + eta.v1~
    c(0.000486,
      -0.000534, 0.0139,
      -0.00258, 0.0225, 0.269,
      0.000845, -0.00527, -0.00809, 0.00751,
      0.0263, -0.418, -0.479, 0.32, 0.208,
      0.00039, -0.00545, -0.00607, 0.00611, 0.335, 0.0108,
      -0.000606, 0.000954, 0.0151, 0.00167, 0.104, 0.00836, 0.0309,
      -0.000145, -0.000984, -0.0166, -0.000275, 0.0548, 0.000295, 0.0113, 0.0284,
      -0.0000000274, 0.00000000391, 0.00000127, -0.000000134, 0.000000296, -0.000000204, 0.0000000143, 0.000000448, 0.000000000128,
      0.000105, -0.000714, -0.0137, -0.000241, 0.00985, -0.000194, -0.00118, 0.00253, -0.0000000584, 0.00511,
      0.000228, 0.000948, 0.00705, -0.0013, -0.0379, -0.000816, -0.00151, -0.00753, -0.000000197, -0.00132, 0.024,
      0.000958, -0.00452, -0.0108, 0.00556, 0.00729, 0.0025, -0.00871, 0.00649, -0.00000207, 0.00492, -0.00363, 0.356,
      0.000396, -0.00313, -0.00614, 0.00148, 0.081, -0.000111, -0.00219, 0.000206, 0.000000026, 0.000287, -0.0000323, 0.00203, 0.00184,
      -0.000215, 0.000461, 0.00443, -0.000661, -0.0237, -0.000494, 0.000357, 0.000208, 0.000000051, -0.0000936, 0.000192, -0.00116, -0.000179, 0.000026
      )
)


#create patients data set
library(tmvtnorm)
#> Loading required package: mvtnorm
#> Loading required package: Matrix
#> Loading required package: stats4
#> Loading required package: gmm
#> Loading required package: sandwich
set.seed(1234)

## 1. define variables
cov_names <- c("BW",                  # Body weight (kg)
               "BMI",                 # Body mass index (kg/m2)
               "FB",                  # 24h fluid balance (L)
               "EGFR",                # eGFR (mL/min/1.73m2)
               "SODI",                # Serum sodium (mmol/L)
               "TOPR")                # Serum total protein (g/L)

## 2. median, max and min
mu     <- c(70.0,   24.7,   2.8,   53.0,   138.4,  66.0)
lower  <- c(36.0,   13.7,  -1.5,    6.0,   121.9,  26.0)
upper  <- c(115.0,  39.6,  13.9,  200.0,   155.9,  90.0)

## 3. covariates of matrix
Sigma <- matrix(
  c(256.9,   69.4,   4.3,  -153.7,   -9.3,   24.8,
    69.4,   25.1,   1.1,   -60.7,   -1.3,   11.1,
    4.3,    1.1,   5.1,   -18.8,    0.9,    0.9,
    -153.7,  -60.7, -18.8,  1375.8,  -31.4,  -98.4,
    -9.3,   -1.3,   0.9,   -31.4,   40.5,   -1.4,
    24.8,   11.1,   0.9,   -98.4,   -1.4,  123.2),
  nrow = 6, byrow = TRUE
)
colnames(Sigma) <- cov_names
rownames(Sigma) <- cov_names

## 4.extract 200 patients
n_virt <- 200

X_trans <- rtmvnorm(
  n      = n_virt,
  mean   = mu,
  sigma  = Sigma,
  lower  = lower,
  upper  = upper
)
colnames(X_trans) <- cov_names

## 5. create virtual patients data frame
virtual_patients <- as.data.frame(X_trans)
virtual_patients$id <- 1:n_virt
virtual_patients <- virtual_patients[, c("id", cov_names)]

#dose and event design
s <- seq(0, 24, by = 8)

#create event table
e <- et() %>%
  et(id=virtual_patients$id, amt=15*virtual_patients$BW, duration=0.5, rate=(15*virtual_patients$BW)/0.5, time=0) %>%
  et(s) %>%
  as.data.frame %>%
  merge(virtual_patients, by="id") %>%
  as_tibble

thetaMat <- as.matrix(nearPD(thetaMat)$mat)

#create simulation
sim <- rxSolve(
  rx1, theta, e,
  thetaMat = thetaMat,
  thetaLower=0,
  omega = c("eta.cl", "eta.v1"),
  omegaXform = "variance",
  sigma = sigma,
  sigmaXform = "identity",
  nStud = 100,
  dfSub=199,
  dfObs=599
)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:01

print(sim)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 20,000 × 16
#>    sim.id id     tvcl tvcl_egfr tvcl_bmi  eta.cl  tvv1 tvv1_bmi tvv1_topr
#>     <int> <fct> <dbl>     <dbl>    <dbl>   <dbl> <dbl>    <dbl>     <dbl>
#>  1      1 1      2.18      1.53     1.00 -1.94    16.6    0.576     0.678
#>  2      1 2      2.18      1.53     1.00  0.218   16.6    0.576     0.678
#>  3      1 3      2.18      1.53     1.00 -0.497   16.6    0.576     0.678
#>  4      1 4      2.18      1.53     1.00  0.0484  16.6    0.576     0.678
#>  5      1 5      2.18      1.53     1.00 -0.954   16.6    0.576     0.678
#>  6      1 6      2.18      1.53     1.00  0.458   16.6    0.576     0.678
#>  7      1 7      2.18      1.53     1.00 -0.0513  16.6    0.576     0.678
#>  8      1 8      2.18      1.53     1.00  1.14    16.6    0.576     0.678
#>  9      1 9      2.18      1.53     1.00 -0.759   16.6    0.576     0.678
#> 10      1 10     2.18      1.53     1.00  0.416   16.6    0.576     0.678
#> # ℹ 19,990 more rows
#> # ℹ 7 more variables: tvv1_sodi <dbl>, tvv1_fb <dbl>, tvv1_flg <dbl>,
#> #   eta.v1 <dbl>, tvq <dbl>, tvv2 <dbl>, prop.sd <dbl>
#> ── Initial Conditions ($inits): ──
#> centr  peri 
#>     0     0 
#> 
#> Simulation with uncertainty in:
#> • parameters ($thetaMat for changes)
#> • omega matrix ($omegaList)
#> • sigma matrix ($sigmaList)
#> 
#> ── First part of data (object): ──
#> # A tibble: 80,000 × 18
#>   sim.id    id  time  TVCL    CL  TVV1    V1     Q    V2 ipred    Dv centr  peri
#>    <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1      1     1     0 3.48  0.502  2.20  3.21  1.15  12.7   0     0     0      0 
#> 2      1     1     8 3.48  0.502  2.20  3.21  1.15  12.7  34.6  27.9 111.   539.
#> 3      1     1    16 3.48  0.502  2.20  3.21  1.15  12.7  25.9  34.8  83.3  449.
#> 4      1     1    24 3.48  0.502  2.20  3.21  1.15  12.7  21.3  19.2  68.4  370.
#> 5      1     2     0 0.442 0.549 23.9  11.2   1.15  12.7   0     0     0      0 
#> 6      1     2     8 0.442 0.549 23.9  11.2   1.15  12.7  42.8  45.8 479.   398.
#> # ℹ 79,994 more rows
#> # ℹ 5 more variables: BMI <dbl>, FB <dbl>, EGFR <dbl>, SODI <dbl>, TOPR <dbl>

^{Created on 2025-12-17 with reprex v2.1.1}

0 replies

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

There is error in simulation with covariance matrix, I don't know what's the problem is. #972

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 2 comments

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

There is error in simulation with covariance matrix, I don't know what's the problem is. #972

Uh oh!

Uh oh!

WXX155367 Dec 17, 2025

Replies: 2 comments

Uh oh!

Uh oh!

WXX155367 Dec 17, 2025 Author

Uh oh!

mattfidler Dec 17, 2025 Maintainer

WXX155367
Dec 17, 2025

WXX155367
Dec 17, 2025
Author

mattfidler
Dec 17, 2025
Maintainer