Feature Request: nlmixr2.formula() function

[billdenney]

As initially implemented in #5, this is an issue for discussion of the feature of nlmixr2.formula().

The overall interface is below. (This isn't going to be precise at the moment as it's for discussion. The final interface will match the generic of nlmixr2().

nlmixr2.formula(object, data, ..., start, param, residualModel=~add(add_err))

• object will be the formula, and upon reconsideration, I think it should be a simple formula without the double tilde notation originally proposed in Add nlmixr_formula() (testable work in progress) #5 (more on that later).
• data used for evaluation of factor levels, possibly modified, and sent to nlmixr2.function().
• ... passed to nlmixr2.function()
• start initial parameter estimates as either a named vector or as a named list. If you have a single estimate per fixed effect, then the named vector and named list are equivalent. If you want to start with different fixed effects per factor level in your data, you need to use a list.
• param is used to define random effects and factor-modeled fixed effects (the information after the second tilde in Add nlmixr_formula() (testable work in progress) #5 moved here).
• residualModel has the residual error model which is sent to model() as the right-hand-side of the prediction.

The most complex element is that param now defines both fixed and random effects based on the way it will be used. This is the way that brms defines its fixed and random effects (see example 3 in https://cran.r-project.org/web/packages/brms/vignettes/brms_multilevel.pdf). (With minimal modification, we could even make this function take in a brms::bf() object as the object argument where that would encompass the object and param arguments.)

An example call would be:

  nlmixr2.formula(
    height ~ Asym+(R0-Asym)*exp(-exp(lrc)*age),
    data = Loblolly,
    start = list(Asym = 103, R0 = -8.5, lrc = -3.3, add_err=1),
    param =
      list(
        Asym~species+(1|Seed),
        lrc~1|Seed
      ),
    est=NULL
  )

(This example doesn't actually work since "species" isn't in the Loblolly data. It is just to show an example.)

This would specify a model where Asym has a separate fixed effect estimated per level of the Loblolly$species factor. There are random effects for Asym and lrc on the Seed identifier.

This is very close to what was just discussed in the team meeting earlier. The only notable difference is moving the random effect specification into param. The benefits of moving it to param are:

• There is a consistent interface to creating parameters rather than some parameters created after an unusual second tilde and others created within param.
• It simplifies the way that parameters can be considered all in one place (so, there should be fewer user errors).
• Generating appropriate mu-referenced covariates should be easier. (i.e. the code should be easier to maintain.)

[billdenney]

Notes from team meeting:

• Confirm that there is only one random effect level.
• Check to make sure that generated names are not in the model formula.
• Allow factor values to be named
• Make vignette
  • using "rxSolve" with an Emax model (with random effect)
  • without random effects
  • put a random effect in later

model({
Asym <- Asym.oak + Asym.elm*(species == "elm") + Asym_re
height <- Asym+(R0-Asym)*exp(-exp(lrc)*age)
height ~ add(add_err)
})

nlmixr2.formula(
    height ~ Asym+(R0-Asym)*exp(-exp(lrc)*age),
    data = Loblolly,
    start = list(Asym = c(oak=103, elm=1), R0 = -8.5, lrc = -3.3, add_err=1),
    param =
      list(
        Asym~species+(1|Seed),
        lrc~1|Seed
      ),
    est=NULL
  )