nlmixrdevelopment / nlmixr

nlmixr: an R package for population PKPD modeling
https://nlmixrdevelopment.github.io/nlmixr/
GNU General Public License v2.0
114 stars 45 forks source link

Cannot update a model with `ini()` #520

Closed billdenney closed 3 years ago

billdenney commented 3 years ago

I'm trying to write a likelihood profiling function, and so, I'm trying to modify the model via ini(), but it gives the following error. I think that I'm doing it right based on the vignette (https://nlmixrdevelopment.github.io/nlmixr/articles/modelPiping.html), but if not, I'm not sure what to do differently:

library(nlmixr)
pheno <- function() {
  ini({
    tcl <- log(0.008) # typical value of clearance
    tv <-  log(0.6)   # typical value of volume
    ## var(eta.cl)
    eta.cl + eta.v ~ c(1, 
                       0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
    # interindividual variability on clearance and volume
    add.err <- 0.1    # residual variability
  })
  model({
    cl <- exp(tcl + eta.cl) # individual value of clearance
    v <- exp(tv + eta.v)    # individual value of volume
    ke <- cl / v            # elimination rate constant
    d/dt(A1) = - ke * A1    # model differential equation
    cp = A1 / v             # concentration in plasma
    cp ~ add(add.err)       # define error model
  })
}

fitted <-
  nlmixr(
    pheno, pheno_sd, "saem",
    control=list(print=0), 
    table=list(cwres=TRUE, npde=TRUE)
  )
#> i parameter labels from comments will be replaced by 'label()'
#> > generate SAEM model
#> v done
#> RxODE 1.0.8.1 using 4 threads (see ?getRxThreads)
#> Calculating covariance matrix
#> Calculating residuals/tables
#> Compiling NPDE model...done
#> done

ini(fitted, tcl <- fixed(-5))
#> Error in fixed(-5): could not find function "fixed"

Created on 2021-04-29 by the reprex package (v2.0.0)

mattfidler commented 3 years ago

Hi @billdenney

This error is because you used the <- instead of the = operator. This is one of the cases where these operators produce different things (since it is within a function).

Here is the modified reprex:

library(nlmixr)

pheno <- function() {
  ini({
    tcl <- log(0.008) # typical value of clearance
    tv <-  log(0.6)   # typical value of volume
    ## var(eta.cl)
    eta.cl + eta.v ~ c(1,
                       0.01, 1) ## cov(eta.cl, eta.v), var(eta.v)
    # interindividual variability on clearance and volume
    add.err <- 0.1    # residual variability
  })
  model({
    cl <- exp(tcl + eta.cl) # individual value of clearance
    v <- exp(tv + eta.v)    # individual value of volume
    ke <- cl / v            # elimination rate constant
    d/dt(A1) = - ke * A1    # model differential equation
    cp = A1 / v             # concentration in plasma
    cp ~ add(add.err)       # define error model
  })
}

fitted <-
  nlmixr(
    pheno, pheno_sd, "saem",
    control=list(print=0),
    table=list(cwres=TRUE, npde=TRUE)
  )
#> ℹ parameter labels from comments will be replaced by 'label()'
#> RxODE 1.0.8.1 using 4 threads (see ?getRxThreads)
#> Calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> Calculating residuals/tables
#> Compiling NPDE model...done
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> done

ini(fitted, tcl=fixed(-5))
#> ▂▂ RxODE-based ODE model ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ 
#> ── Initialization: ───────────────────────────────────────────────────────────── 
#> Fixed Effects ($theta): 
#>        tcl         tv 
#> -5.0000000  0.3412029 
#> 
#> Omega ($omega): 
#>           eta.cl     eta.v
#> eta.cl 0.2269720 0.1843463
#> eta.v  0.1843463 0.1626408
#> ── μ-referencing ($muRefTable): ──────────────────────────────────────────────── 
#>                              ┌─────────┬─────────┐
#>                              │ theta   │ eta     │
#>                              ├─────────┼─────────┤
#>                              │ tcl     │ eta.cl  │
#>                              ├─────────┼─────────┤
#>                              │ tv      │ eta.v   │
#>                              └─────────┴─────────┘
#> 
#> ── Model: ────────────────────────────────────────────────────────────────────── 
#>     cl <- exp(tcl + eta.cl) # individual value of clearance
#>     v <- exp(tv + eta.v)    # individual value of volume
#>     ke <- cl / v            # elimination rate constant
#>     d/dt(A1) = - ke * A1    # model differential equation
#>     cp = A1 / v             # concentration in plasma
#>     cp ~ add(add.err)       # define error model 
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂

Created on 2021-05-01 by the reprex package (v2.0.0)

mattfidler commented 3 years ago

Thank you for your patience in this response; I'm a bit swamped with other projects at the moment.

billdenney commented 3 years ago

Thanks!