Simulating models with uncertainty that have covariates

[cfosser]

Reminder from ACoP9: Please add functionality of simulating models that have covariates added. Specific case is in the hands on exercise called 'Case Theophylline' - after Step 4 - cannot simulate from this model fit.

[mattfidler]

See https://github.com/nlmixrdevelopment/RxODE/issues/26

[billdenney]

I was about to ask a similar question about how to simulate. It would be helpful to add it as a section to the website.

Simulations from nlme have an interesting issue where you will need to ensure that the subject identifiers differ so that you can ensure that you're getting a true random simulation (perhaps best would be a warning message when subject numbers overlap with the original estimation subject numbers).

[mattfidler]

Can you give an example of what you mean Bill

[billdenney]

@mattfidler, Apologies, I was remembering predictions from lme4 which can have issues when predicting into a mix of new and old levels. As I recall, there is no (good) way to predict nlme models with new data without writing your own prediction function.

Revised below to show one more issue: new levels just get the population prediction.

In the case below (with lme4, so not applicable...), subjects 308 and 309 maintain their random effects from the underlying model while subjects 400 and 401 gets no random effect. If subject identifiers are integers starting at 1 or any other scenario where the subject identifier may overlap, you may not get the randomness you expect from the new subject numbering. I've been bitten with this a few times in the past (and I now have a best practice where my simulated subject identifiers add 1e6 so that they are very unlikely to overlap with subjects estimated random effects from the underlying model).

The specific note below is that the first two predicted values are identical when I may have intended them to differ.

library(lme4)
#> Loading required package: Matrix
fm1 <- lmer(Reaction ~ Days + (Days | Subject), sleepstudy)
unique(sleepstudy$Subject)
#>  [1] 308 309 310 330 331 332 333 334 335 337 349 350 351 352 369 370 371
#> [18] 372
#> 18 Levels: 308 309 310 330 331 332 333 334 335 337 349 350 351 352 ... 372
predict(fm1, newdata=data.frame(Subject=factor(c(308, 309)), Days=3))
#>        1        2 
#> 312.6624 216.5493
predict(fm1, newdata=data.frame(Subject=factor(c(308, 309, 400)), Days=3))
#> Error in levelfun(r, n, allow.new.levels = allow.new.levels): new levels detected in newdata
predict(fm1, newdata=data.frame(Subject=factor(c(308, 309, 400, 401)), Days=3), allow.new.levels=TRUE)
#>        1        2        3        4 
#> 312.6624 216.5493 282.8070 282.8070

Created on 2018-11-13 by the reprex package (v0.2.0).

[mattfidler]

Thank you Bill for clarifying.

Simulation needs to be documented better. There is some documentation in the PAGE & ACOP presentations.

https://github.com/nlmixrdevelopment/ACoP9-2018/blob/master/nlmixr-tutorial_final.pdf

Starting on Slide 27 is the most recent discussion about simulating with nlmixr models. Using nlme with the nlmixr functions can be used to simulate based on the model that nlme provided. The same interface is used for all other estimation methods right now.

We are thinking of expanding this to add covariates and other features.

[mattfidler]

The lme4 methods are not yet supported with the UI yet. I think they use the same Bates/Pinhero methodology. However, for some reason, if I recall correctly, it couldn't fit the theo_sd dataset.

@wwang-at-github do you know if this was ever fixed in lme4?

[wwang-at-github]

maybe not. https://github.com/lme4/lme4/issues/363

On Tue, Nov 13, 2018 at 12:15 PM Matthew Fidler notifications@github.com wrote:

The lme4 methods are not yet supported with the UI yet. I think they use the same Bates/Pinhero methodology. However, for some reason, if I recall correctly, it couldn't fit the theo_sd dataset.

@wwang-at-github https://github.com/wwang-at-github do you know if this was ever fixed in lme4?

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/nlmixrdevelopment/nlmixr/issues/84#issuecomment-438355003, or mute the thread https://github.com/notifications/unsubscribe-auth/AMdZdqmeno49SXhcgbchnj2DudW5lb7Pks5uuv5IgaJpZM4XYCWr .

[Potamophylax]

Hi, @mattfidler ! Thank you dear colleagues for this useful discussion!

I have a small question to which I don't have a direct answer.

What should I do, if I want to simulate the trials, in general, corresponding to the project received in nlmixr. The only difference is that, I have covariates. And I want them to accept fixed values. For example, all patients weighing 50 kg.

Please note that if I substitute the value of the desired body weight in the "parameters", the function "simulate" will require me to enter all the other parameters values. But I don't want to take other options. I want to take those values ​​that were obtained using SAEM in nlmixr in the frame of used project.

res <- readRDS(file = "run1.res.rds")
ev <- eventTable(amount.units="mg", time.units="hours") %>%
  add.dosing(dose=240, nbr.doses=5, dosing.interval = 336, dosing.to=1) %>%
  add.sampling(seq(0,1800,length.out=200))

sim <- simulate(res, events = ev, nSub = 100, nStud = 250,
                params=c(BW=70, IGG=1000, iggeff=0.632, bweff=2.229, tv3=1.093, tv2=1.048, tcl=-4.466, tq=-3.389) )

I would to use the corrected values ​​of only 2 covariates in the problem (BW and IGG), and not all the parameters of the model.

[mattfidler]

Currently what you ask for is not possible, but it should be something that could be added.

[Potamophylax]

Thanks, Matt!

I am very glad to hear this: these innovations would make the work with your package more convenient!

[mattfidler]

With the new RxODE, you can simulate individual covariates from a nlmixr model with model piping;

An example of model piping without covariates is found here.

However, it would be easy enough to specify individual non-time varying covariates by using the iCov in the new RxODE rxSolve:

fit %>% 
    et(amt=300,ii=12,until=48) %>% # Now use amt=300 BID for 48 hrs
    rxSolve(nSub=100, iCov=data.frame(WT=rnorm(100,70,3)) %>% # Simulate with 100 subjects
    confint() %>% # Create intervals (not exactly confidence intervals)
    plot() # Plot the intervals

Although I think predict would also make sense, but it hasn't been implemented yet. Nor has it been implemented in broom's augment.