addCov works when simulating from IIV but not parameters dataframe

[john-harrold]

In the example below I'm simulating first using IIV. I can see the covariate (WT) when I run this simulation. Then I use that to create a table of parameter values and the covariate. When I simulate from that table I cannot see the covariate in the output. I think it should be there (based on what I'm doing at least) right?

library(dplyr)
library(nlmixr2lib)
library(rxode2)

mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)*(WT/70)^.75
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)*(WT/70)^.75

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))

iCov = data.frame(id=c(1:3), WT=70*exp(rnorm(mean=0, n=3, sd=.1)))
sim = rxSolve(mymod, events = et(c(0,1), id=c(1:3)), iCov = iCov, nSub=3, addCov=TRUE)

# Here the covariate is present
head(as.data.frame(sim)$WT)

params = sim$params
params$id = as.numeric(params$id)

subjects = dplyr::left_join(params, iCov, by=c("id"))

sim = rxSolve(mymod, params=subjects, events = et(c(0,1), id=c(1:3)), addCov=TRUE)
# Here the covariate is not present:
head(as.data.frame(sim)$WT)

[mattfidler]

Hi John,

Since you add the covariant WT into the params object it becomes an input parameter instead of a covariate.

So using your reprex above with slight modifications we see it is still in the output object:

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(nlmixr2lib)
library(rxode2)
#> rxode2 2.1.1 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)*(WT/70)^.75
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)*(WT/70)^.75

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ promote `etalka` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalka` to `0.1`
#> ℹ promote `etalcl` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalcl` to `0.1`
#> ℹ promote `etalq` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalq` to `0.1`
#> ℹ promote `etalvp` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvp` to `0.1`
#> ℹ promote `etalvc` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvc` to `0.1`

iCov = data.frame(id=c(1:3), WT=70*exp(rnorm(mean=0, n=3, sd=.1)))
sim = rxSolve(mymod, events = et(c(0,1), id=c(1:3)), iCov = iCov, nSub=3, addCov=TRUE)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments

# Here the covariate is present
head(as.data.frame(sim)$WT)
#> [1] 67.44844 67.44844 72.86528 72.86528 57.64116 57.64116

params = sim$params
params$id = as.numeric(params$id)

subjects = dplyr::left_join(params, iCov, by=c("id"))

sim = rxSolve(mymod, params=subjects, events = et(c(0,1), id=c(1:3)), addCov=TRUE)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> Warning: multi-subject simulation without without 'omega'

# you can see the weight in the `params` object:
print(sim)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#> # A tibble: 3 × 13
#>   id      lka   lcl   lvc   lvp    lq propSd  etalka etalcl  etalq  etalvp
#>   <fct> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl>  <dbl>  <dbl>   <dbl>
#> 1 1      0.45     1     3     5   0.1    0.5 -0.504   0.478  0.214 -0.0108
#> 2 2      0.45     1     3     5   0.1    0.5  0.0665 -0.427  0.522  0.655 
#> 3 3      0.45     1     3     5   0.1    0.5  0.333  -0.423 -0.199  0.224 
#> # ℹ 2 more variables: etalvc <dbl>, WT <dbl>
#> ── Initial Conditions ($inits): ──
#>       depot     central peripheral1 
#>           0           0           0 
#> ── First part of data (object): ──
#> # A tibble: 6 × 16
#>      id  time    ka    cl    vc    vp     q    kel    k12     k21    Cc ipredSim
#>   <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>   <dbl> <dbl>    <dbl>
#> 1     1     0 0.948  4.26  21.2  147. 1.33  0.201  0.0627 0.00907     0        0
#> 2     1     1 0.948  4.26  21.2  147. 1.33  0.201  0.0627 0.00907     0        0
#> 3     2     0 1.68   1.83  11.1  286. 1.92  0.165  0.173  0.00672     0        0
#> 4     2     1 1.68   1.83  11.1  286. 1.92  0.165  0.173  0.00672     0        0
#> 5     3     0 2.19   1.54  29.9  186. 0.783 0.0514 0.0261 0.00421     0        0
#> 6     3     1 2.19   1.54  29.9  186. 0.783 0.0514 0.0261 0.00421     0        0
#> # ℹ 4 more variables: sim <dbl>, depot <dbl>, central <dbl>, peripheral1 <dbl>

# Verified here:
print(head(sim$params))
#>   id  lka lcl lvc lvp  lq propSd     etalka     etalcl      etalq      etalvp
#> 1  1 0.45   1   3   5 0.1    0.5 -0.5038562  0.4779739  0.2143264 -0.01075541
#> 2  2 0.45   1   3   5 0.1    0.5  0.0664669 -0.4266709  0.5221115  0.65450235
#> 3  3 0.45   1   3   5 0.1    0.5  0.3330779 -0.4232901 -0.1994148  0.22449166
#>        etalvc       WT
#> 1  0.05634376 67.44844
#> 2 -0.59450683 72.86528
#> 3  0.39912881 57.64116

^{Created on 2023-12-26 with reprex v2.0.2}

[mattfidler]

I suppose it could be retained since the UI declares WT to be a covariate, though the low level rxode2 models have never made this disinction. keep, drop and addCov always related to the input dataset instead of the input parameters.

To do this, a new feature #638 would need to be integrated first.

[john-harrold]

I think I can work with this. The only downside here is if the covariate changes with time. I don't know how to get the interpolated values out of the simulation.

[mattfidler]

I think I can work with this. The only downside here is if the covariate changes with time. I don't know how to get the interpolated values out of the simulation.

Both iCov and parameter-based simulations enforce one covariate per id (meaning time varying covariates don't exist).

However, if you are simulating from a dataset with time-varying covariates, the time varying covariates have to be in the simulation dataset. The simulation dataset can simply be a NONMEM /rxode2 input dataset. When simulating with resampling (which is highly suggested). you can use resample=TRUE and the corresponding covariate interpolation/simulation for time varying covariates will be in the output dataset.

It is a bit cumbersome to be in 2 different places; Using the iCov option merges the individual covariates with the input dataset and makes a single place to check for covariates (ie the output).

Also using returnType="data.frame", returnType="tibble" or `returnType="data.table" will not drop the individual covariate. (Note that these sorts of solves are slightly faster than the default return type).

[mattfidler]

Hopefully this helps answer you question.

[john-harrold]

However, if you are simulating from a dataset with time-varying covariates, the time varying covariates have to be in the simulation dataset. The simulation dataset can simply be a NONMEM /rxode2 input dataset. When simulating with resampling (which is highly suggested). you can use resample=TRUE and the corresponding covariate interpolation/simulation for time varying covariates will be in the output dataset.

Can you point me to an example of this where I can see what the structure looks like? Sorry if I'm kind of slow here, I've just never used NONMEM for simulation because I just don't enjoy pain that much :)

[mattfidler]

Here is am example:

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(nlmixr2lib)
library(rxode2)
#> rxode2 2.1.1 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

set.seed(42)

mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)*(WT/70)^.75
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)*(WT/70)^.75

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ promote `etalka` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalka` to `0.1`
#> ℹ promote `etalcl` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalcl` to `0.1`
#> ℹ promote `etalq` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalq` to `0.1`
#> ℹ promote `etalvp` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvp` to `0.1`
#> ℹ promote `etalvc` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvc` to `0.1`

iCov = data.frame(id=c(1:3), WT0=70*exp(rnorm(mean=0, n=3, sd=.1)))

events <- et(seq(0.5,12,length.out=10)) |> et(amt=10) |> et(id=1:3) |> as.data.frame()

events <- dplyr::inner_join(events, iCov) |>
  dplyr::mutate(WT=WT0 + 0.1 * time)
#> Joining with `by = join_by(id)`

sim = rxSolve(mymod, events = events)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
head(sim)
#>   id     time       ka       cl       vc       vp        q       kel        k12
#> 1  1 0.500000 2.263795 3.175871 28.68113 121.0086 1.633486 0.1107303 0.05695332
#> 2  1 1.777778 2.263795 3.179658 28.68113 121.0086 1.635434 0.1108624 0.05702125
#> 3  1 3.055556 2.263795 3.183445 28.68113 121.0086 1.637381 0.1109944 0.05708915
#> 4  1 4.333333 2.263795 3.187230 28.68113 121.0086 1.639328 0.1111264 0.05715702
#> 5  1 5.611111 2.263795 3.191013 28.68113 121.0086 1.641274 0.1112583 0.05722487
#> 6  1 6.888889 2.263795 3.194795 28.68113 121.0086 1.643219 0.1113901 0.05729268
#>          k21        Cc  ipredSim       sim        depot  central peripheral1
#> 1 0.01349892 0.2248754 0.2248754 0.1287080 3.224211e+00 6.449680   0.1105820
#> 2 0.01351502 0.2729823 0.2729823 0.4520469 1.787190e-01 7.829440   0.6719465
#> 3 0.01353111 0.2258465 0.2258465 0.3684850 9.906455e-03 6.477534   1.1778044
#> 4 0.01354720 0.1832342 0.1832342 0.1440957 5.491165e-04 5.255364   1.5802540
#> 5 0.01356328 0.1487477 0.1487477 0.1255200 3.043942e-05 4.266251   1.8964349
#> 6 0.01357935 0.1210405 0.1210405 0.1286250 1.687647e-06 3.471580   2.1431771
#>         WT
#> 1 80.33566
#> 2 80.46344
#> 3 80.59122
#> 4 80.71900
#> 5 80.84678
#> 6 80.97455

plot(sim, Cc)

^{Created on 2023-12-26 with reprex v2.0.2}

You simply add it to the input dataset or the event dataset. In this case the WT is a time changing covariate.

[mattfidler]

And here is an example where you resample the covariates:

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(nlmixr2lib)
library(rxode2)
#> rxode2 2.1.1 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

set.seed(42)

mymod = function() {
  description <- "Two compartment PK model with linear clearance using differential equations"
  ini({
    lka <- 0.45 ; label("Absorption rate (Ka)")
    lcl <- 1 ; label("Clearance (CL)")
    lvc <- 3 ; label("Central volume of distribution (V)")
    lvp <- 5 ; label("Peripheral volume of distribution (Vp)")
    lq  <- 0.1 ; label("Intercompartmental clearance (Q)")
    propSd <- 0.5 ; label("Proportional residual error (fraction)")
  })
  model({
    ka <- exp(lka)
    cl <- exp(lcl)*(WT/70)^.75
    vc <- exp(lvc)
    vp <- exp(lvp)
    q  <- exp(lq)*(WT/70)^.75

    kel <- cl/vc
    k12 <- q/vc
    k21 <- q/vp

    d/dt(depot) <- -ka*depot
    d/dt(central) <-  ka*depot - kel*central - k12*central + k21*peripheral1
    d/dt(peripheral1) <- k12*central - k21*peripheral1
    Cc <- central / vc

    Cc ~ prop(propSd)
  })
}

mymod = addEta(mymod, c("lka", 'lcl', 'lq', 'lvp', 'lvc'))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ promote `etalka` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalka` to `0.1`
#> ℹ promote `etalcl` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalcl` to `0.1`
#> ℹ promote `etalq` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalq` to `0.1`
#> ℹ promote `etalvp` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvp` to `0.1`
#> ℹ promote `etalvc` to between subject variability with initial estimate 0.1
#> ℹ change initial estimate of `etalvc` to `0.1`

iCov = data.frame(id=c(1:3), WT0=70*exp(rnorm(mean=0, n=3, sd=.1)))

events <- et(seq(0.5,12,length.out=10)) |> et(amt=10) |> et(id=1:3) |> as.data.frame()

events <- dplyr::inner_join(events, iCov) |>
  dplyr::mutate(WT=WT0 + 0.1 * time)
#> Joining with `by = join_by(id)`

sim = rxSolve(mymod, events = events,
              nStud=4, resample=TRUE)
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
head(sim)
#>   sim.id id     time       ka       cl       vc       vp        q        kel
#> 1      1  1 0.500000 2.263795 2.746994 28.68113 121.0086 1.412896 0.09577705
#> 2      1  1 1.777778 2.263795 2.750969 28.68113 121.0086 1.414941 0.09591565
#> 3      1  1 3.055556 2.263795 2.754943 28.68113 121.0086 1.416985 0.09605419
#> 4      1  1 4.333333 2.263795 2.758914 28.68113 121.0086 1.419027 0.09619266
#> 5      1  1 5.611111 2.263795 2.762884 28.68113 121.0086 1.421069 0.09633106
#> 6      1  1 6.888889 2.263795 2.766851 28.68113 121.0086 1.423110 0.09646940
#>          k12        k21        Cc  ipredSim       sim        depot  central
#> 1 0.04926222 0.01167599 0.2263708 0.2263708 0.1295639 3.224211e+00 6.492571
#> 2 0.04933351 0.01169289 0.2814011 0.2814011 0.4659881 1.787186e-01 8.070902
#> 3 0.04940476 0.01170978 0.2393803 0.2393803 0.3905663 9.906402e-03 6.865699
#> 4 0.04947598 0.01172666 0.1997851 0.1997851 0.1571113 5.491133e-04 5.730064
#> 5 0.04954717 0.01174353 0.1667754 0.1667754 0.1407326 3.043808e-05 4.783307
#> 6 0.04961832 0.01176040 0.1394626 0.1394626 0.1482013 1.687227e-06 3.999944
#>   peripheral1       WT
#> 1  0.09605911 66.20665
#> 2  0.59104984 66.33443
#> 3  1.04958458 66.46221
#> 4  1.42622621 66.58998
#> 5  1.73236286 66.71776
#> 6  1.98011910 66.84554

plot(sim, Cc)

^{Created on 2023-12-26 with reprex v2.0.2}

[john-harrold]

What I'm trying to understand is the scenario where you have a set of parameters for subjects with observations of a time-varying covariate. I think should be able to pass the parameters to this call:

sim = rxSolve(mymod, events = events, nStud=4, resample=TRUE)

What do I do when you have observations in events but do not have covariant values. Like if you have weight measurements at each cycle but you want to simulate at times during the cycle? Can I just set the covariate values at those times to NA?

[mattfidler]

That would work. NA values are interpolated using whatever interpolation method you are using throughout (ie locf nocb) the output should still show the interpolated values (not NA)