Give a better error with separation strategy that is missing the `dfObs` and `dfSub`

[mattfidler]

    Thanks... this brought me a couple of step forward :) Unfortunately, on the model I am working on (which, regrettably, I cannot share it due to company policies) I get a warning `chol(): given matrix is not symmetric` when I simulate with the choleski matrix. This means that there are NaNs in the simulated `omegaList`. The model random effects are expressed as: `cl = exp(tcl + eta.cl) ` (and I am reading the cov file from NONMEM results). Is the correct `omegaXform` argument  still `variance` or should it be `log` or something else?  I do not get the same issue if I set `omegaXform = identity`, but the results do not look reasonable. I will see if I can recreate this in a minimal working example, but is there a way that I can further debug/dig into this?  I completely understand that it is difficult to say without the model and the parameter values...

I have also found another nice-to-fix. I was being stupid and I was attempting to simulate uncertainty in the population parameters forgetting the dfSub and dfObs arguments in rxSolve, see the minmal example below. The error message (Error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits, : matrix must be a square matrix) is misleading and it would be nice to have an error message similar to dfSub and dfObs must be specified.

rm(list = ls())

library(tidyverse)
library(nlmixr2)
#> Loading required package: nlmixr2data

rx1 <- RxODE({
  cl <- tcl*(1+crcl.cl*(CLCR-65)) * exp(eta.v)
  v <- tv * WT * exp(eta.v)
  ka <- tka * exp(eta.ka)
  ipred <- linCmt()
  obs <- ipred * (1 + prop.sd) + add.sd 
})

theta <- c(tcl=2.63E+01, tv=1.35E+00, tka=4.20E+00, tlag=2.08E-01,
           prop.sd=2.05E-01, add.sd=1.06E-02, crcl.cl=7.17E-03,
           ## Note that since we are using the separation strategy the ETA variances are here too
           eta.cl=7.30E-02,  eta.v=3.80E-02, eta.ka=1.91E+00)

thetaMat <- lotri(
  tcl + tv + tka + tlag + prop.sd + add.sd + crcl.cl + eta.cl + eta.v + eta.ka ~
    c(7.95E-01,
      ##      2.05E-02, 1.92E-03, --> Here I am assuming that the tv was fixed during estimation
      ##    so that nonmem cov output is a low triangular matrix with a zero row
      0, 0,
      7.22E-02, -8.30E-03, 6.55E-01,
      -3.45E-03, -6.42E-05, 3.22E-03, 2.47E-04,
      8.71E-04, 2.53E-04, -4.71E-03, -5.79E-05, 5.04E-04,
      6.30E-04, -3.17E-06, -6.52E-04, -1.53E-05, -3.14E-05, 1.34E-05,
      -3.30E-04, 5.46E-06, -3.15E-04, 2.46E-06, 3.15E-06, -1.58E-06, 2.88E-06,
      -1.29E-03, -7.97E-05, 1.68E-03, -2.75E-05, -8.26E-05, 1.13E-05, -1.66E-06, 1.58E-04,
      -1.23E-03, -1.27E-05, -1.33E-03, -1.47E-05, -1.03E-04, 1.02E-05, 1.67E-06, 6.68E-05, 1.56E-04,
      7.69E-02, -7.23E-03, 3.74E-01, 1.79E-03, -2.85E-03, 1.18E-05, -2.54E-04, 1.61E-03, -9.03E-04, 3.12E-01))

thetaMat1 <- thetaMat

w <- which(dimnames(thetaMat)[[1]] == "tv")
thetaMat1 <- thetaMat[-w, -w]

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

thetaMatChol <- chol(thetaMat1)

evw <- et(amount.units="mg", time.units="hours") %>%
  et(amt=100) %>%
  ## For this problem we will simulate with sampling windows
  et(list(c(0, 0.5),
          c(0.5, 1),
          c(1, 3),
          c(3, 6),
          c(6, 12))) %>%
  et(id=1:1000)

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

sim  <- rxSolve(rx1, theta, evw,  nSub=100, nStud=10,
                thetaMat=thetaMatChol,
                thetaIsChol=TRUE,
                ## Match boundaries of problem
                thetaLower=0, 
                sigma=c("prop.sd", "add.sd"), ## Sigmas are standard deviations
                sigmaXform="identity", # default sigma xform="identity"
                omega=c("eta.cl", "eta.v", "eta.ka"), ## etas are variances
                omegaXform="variance", # default omega xform="variance"
                iCov=data.frame(WT=rnorm(1000, 70, 15), CLCR=rnorm(1000, 65, 25)),
                # ! Stupidly forgetting dfSub and dfObs
                # dfSub=74, dfObs=476
                );
#> Error in rxSolveSEXP(object, .ctl, .nms, .xtra, params, events, inits, : matrix must be a square matrix

^{Created on 2022-10-17 by the reprex package (v2.0.1)}

Session info

``` r sessioninfo::session_info() #> ─ Session info ─────────────────────────────────────────────────────────────── #> setting value #> version R version 4.2.0 (2022-04-22) #> os Fedora 33 (Container Image) #> system x86_64, linux-gnu #> ui X11 #> language en_US:en #> collate en_US.UTF-8 #> ctype en_US.UTF-8 #> tz Europe/Zurich #> date 2022-10-17 #> pandoc 2.17.1.1 @ /usr/lib/rstudio-server/bin/quarto/bin/ (via rmarkdown) #> #> ─ Packages ─────────────────────────────────────────────────────────────────── #> package * version date (UTC) lib source #> assertthat 0.2.1 2019-03-21 [2] CRAN (R 4.2.0) #> backports 1.4.1 2021-12-13 [2] CRAN (R 4.2.0) #> broom 0.8.0 2022-04-13 [2] CRAN (R 4.2.0) #> cachem 1.0.6 2021-08-19 [2] CRAN (R 4.2.0) #> cellranger 1.1.0 2016-07-27 [2] CRAN (R 4.2.0) #> checkmate 2.1.0 2022-04-21 [2] CRAN (R 4.2.0) #> cli 3.3.0 2022-04-25 [2] CRAN (R 4.2.0) #> clipr 0.8.0 2022-02-22 [2] CRAN (R 4.2.0) #> colorspace 2.0-3 2022-02-21 [2] CRAN (R 4.2.0) #> crayon 1.5.1 2022-03-26 [2] CRAN (R 4.2.0) #> data.table 1.14.2 2021-09-27 [2] CRAN (R 4.2.0) #> DBI 1.1.2 2021-12-20 [2] CRAN (R 4.2.0) #> dbplyr 2.2.1 2022-06-27 [1] CRAN (R 4.2.0) #> digest 0.6.29 2021-12-01 [2] CRAN (R 4.2.0) #> dparser 1.3.1-6 2022-10-05 [1] CRAN (R 4.2.0) #> dplyr * 1.0.9 2022-04-28 [2] CRAN (R 4.2.0) #> ellipsis 0.3.2 2021-04-29 [2] CRAN (R 4.2.0) #> evaluate 0.15 2022-02-18 [2] CRAN (R 4.2.0) #> fansi 1.0.3 2022-03-24 [2] CRAN (R 4.2.0) #> fastmap 1.1.0 2021-01-25 [2] CRAN (R 4.2.0) #> forcats * 0.5.1 2021-01-27 [2] CRAN (R 4.2.0) #> fs 1.5.2 2021-12-08 [2] CRAN (R 4.2.0) #> gargle 1.2.0 2021-07-02 [2] CRAN (R 4.2.0) #> generics 0.1.2 2022-01-31 [2] CRAN (R 4.2.0) #> ggplot2 * 3.3.6 2022-05-03 [2] CRAN (R 4.2.0) #> glue 1.6.2 2022-02-24 [2] CRAN (R 4.2.0) #> googledrive 2.0.0 2021-07-08 [2] CRAN (R 4.2.0) #> googlesheets4 1.0.0 2021-07-21 [2] CRAN (R 4.2.0) #> gtable 0.3.0 2019-03-25 [2] CRAN (R 4.2.0) #> haven 2.5.0 2022-04-15 [2] CRAN (R 4.2.0) #> highr 0.9 2021-04-16 [2] CRAN (R 4.2.0) #> hms 1.1.1 2021-09-26 [2] CRAN (R 4.2.0) #> htmltools 0.5.2 2021-08-25 [2] CRAN (R 4.2.0) #> httr 1.4.3 2022-05-04 [2] CRAN (R 4.2.0) #> jsonlite 1.8.0 2022-02-22 [2] CRAN (R 4.2.0) #> knitr 1.39 2022-04-26 [2] CRAN (R 4.2.0) #> lattice 0.20-45 2021-09-22 [2] CRAN (R 4.2.0) #> lbfgsb3c 2020-3.2 2020-03-03 [1] CRAN (R 4.2.0) #> lifecycle 1.0.1 2021-09-24 [2] CRAN (R 4.2.0) #> lotri 0.4.2 2022-06-18 [1] CRAN (R 4.2.0) #> lubridate 1.8.0 2021-10-07 [2] CRAN (R 4.2.0) #> magrittr 2.0.3 2022-03-30 [2] CRAN (R 4.2.0) #> Matrix 1.5-1 2022-09-13 [1] CRAN (R 4.2.0) #> memoise 2.0.1 2021-11-26 [2] CRAN (R 4.2.0) #> modelr 0.1.8 2020-05-19 [2] CRAN (R 4.2.0) #> munsell 0.5.0 2018-06-12 [2] CRAN (R 4.2.0) #> n1qn1 6.0.1-10 2020-11-17 [1] CRAN (R 4.2.0) #> nlme 3.1-157 2022-03-25 [2] CRAN (R 4.2.0) #> nlmixr2 * 2.0.7 2022-06-27 [1] CRAN (R 4.2.0) #> nlmixr2data * 2.0.7 2022-04-22 [1] CRAN (R 4.2.0) #> nlmixr2est 2.0.8 2022-06-22 [1] CRAN (R 4.2.0) #> nlmixr2plot 2.0.6 2022-05-23 [1] CRAN (R 4.2.0) #> pillar 1.7.0 2022-02-01 [2] CRAN (R 4.2.0) #> pkgconfig 2.0.3 2019-09-22 [2] CRAN (R 4.2.0) #> PreciseSums 0.5 2022-04-22 [1] CRAN (R 4.2.0) #> purrr * 0.3.4 2020-04-17 [2] CRAN (R 4.2.0) #> qs 0.25.4 2022-08-09 [1] CRAN (R 4.2.0) #> R6 2.5.1 2021-08-19 [2] CRAN (R 4.2.0) #> RApiSerialize 0.1.2 2022-08-25 [1] CRAN (R 4.2.0) #> Rcpp 1.0.8.3 2022-03-17 [2] CRAN (R 4.2.0) #> RcppParallel 5.1.5 2022-01-05 [2] CRAN (R 4.2.0) #> readr * 2.1.2 2022-01-30 [2] CRAN (R 4.2.0) #> readxl 1.4.0 2022-03-28 [2] CRAN (R 4.2.0) #> reprex 2.0.1 2021-08-05 [2] CRAN (R 4.2.0) #> rlang 1.0.2 2022-03-04 [2] CRAN (R 4.2.0) #> rmarkdown 2.14 2022-04-25 [2] CRAN (R 4.2.0) #> rstudioapi 0.13 2020-11-12 [2] CRAN (R 4.2.0) #> rvest 1.0.2 2021-10-16 [2] CRAN (R 4.2.0) #> rxode2 2.0.7 2022-05-17 [1] CRAN (R 4.2.0) #> scales 1.2.0 2022-04-13 [2] CRAN (R 4.2.0) #> sessioninfo 1.2.2 2021-12-06 [2] CRAN (R 4.2.0) #> stringfish 0.15.7 2022-04-13 [1] CRAN (R 4.2.0) #> stringi 1.7.6 2021-11-29 [2] CRAN (R 4.2.0) #> stringr * 1.4.0 2019-02-10 [2] CRAN (R 4.2.0) #> symengine 0.2.1 2022-04-28 [1] CRAN (R 4.2.0) #> sys 3.4 2020-07-23 [2] CRAN (R 4.2.0) #> tibble * 3.1.7 2022-05-03 [2] CRAN (R 4.2.0) #> tidyr * 1.2.0 2022-02-01 [2] CRAN (R 4.2.0) #> tidyselect 1.1.2 2022-02-21 [2] CRAN (R 4.2.0) #> tidyverse * 1.3.2 2022-07-18 [1] CRAN (R 4.2.0) #> tzdb 0.3.0 2022-03-28 [2] CRAN (R 4.2.0) #> units 0.8-0 2022-02-05 [2] CRAN (R 4.2.0) #> utf8 1.2.2 2021-07-24 [2] CRAN (R 4.2.0) #> vctrs 0.4.1 2022-04-13 [2] CRAN (R 4.2.0) #> vpc 1.2.2 2021-01-11 [1] CRAN (R 4.2.0) #> withr 2.5.0 2022-03-03 [2] CRAN (R 4.2.0) #> xfun 0.30 2022-03-02 [2] CRAN (R 4.2.0) #> xml2 1.3.3 2021-11-30 [2] CRAN (R 4.2.0) #> yaml 2.3.5 2022-02-21 [2] CRAN (R 4.2.0) #> #> [1] /home/brizzif/R/x86_64-pc-linux-gnu-library/4.2.0-foss #> [2] /apps/rocs/2020.08/cascadelake/software/R/4.2.0-foss-2020a/lib64/R/library #> #> ────────────────────────────────────────────────────────────────────────────── ```

Once again, thank you so much!

_Originally posted by @frbrz in https://github.com/nlmixr2/rxode2random/issues/2

[mattfidler]

It will message something like:

  Assertion on 'dfObs' failed: Element 1 is not >= 2.

[frbrz]

Fantastic! Thanks for all the fixes... Looking forward to the next release :)

(and I will try to come back to example with uncertainty in simulation parameters at the end of the week and see if I can reproduce the strange behavior I have for my problem).

[mattfidler]

The next release hit CRAN today for 'rxode2'.

This breaks 'nlmixr2' for now, though.

I have a draft blog post here:

https://deploy-preview-5--nlmixr.netlify.app/blog/2022-10-19-rxode2-2.0.9-release/

[frbrz]

Congrats! The post is very nice. Would it make sense to add to the blog entry a small warning saying that the new rxode2 breaks nlmixr2, and saying that a new nlmixr2 will be available soon (I presume)? And maybe suggest waiting to update rxode2 until nlmixr2 version XYZ is available?

[mattfidler]

It was also available yesterday, but I updated it more to talk about how to proceed in case the binaries are out of sync.

[mattfidler]

Since nlmixr2extra is now failing since it is using an old version of rxode2parse or rxode2, I am resubmitting with version dependencies all the packages needed....

[john-harrold]

Related to nlmixr2extra: will this update have the nlmixr object for testing included :)?

[mattfidler]

It already does, yes. I was about to message you @john-harrold but you beat met to it :smile: