single subject data

[lirui0321]

Hello,

I understand that the same question has been posted previously ( https://github.com/nlmixrdevelopment/nlmixr/issues/177 ), however, the solution provided does not seem to work with nlmixr2. So, what's the right way to fit data with N of 1?

one.cmt <- function() {
    ini({
        tka <- 0.45 # Log Ka
        tcl <- 1 # Log Cl
        tv <- 3.45    # Log V
        eta.ka ~ fix(0)
        add.err <- 0.7
    })
    model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl)
        v <- exp(tv)
        linCmt() ~ add(add.err)
    })
}
single <- nlmixr(one.cmt, theo_sd,"fo",list(print=0))
single <- nlmixr(one.cmt, theo_sd,"bobyqa",list(print=0))
single <- nlmixr(one.cmt, theo_sd,"nls",list(print=0))

I got the following error message

Error in exp(.v) : non-numeric argument to mathematical function
Error: non-numeric argument to mathematical function

[mattfidler]

Hm, I'm not getting this at all:

    library(nlmixr2)
#> Loading required package: nlmixr2data

    one.cmt <- function() {
      ini({
        tka <- 0.45 # Log Ka
        tcl <- 1 # Log Cl
        tv <- 3.45    # Log V
        eta.ka ~ fix(0)
        add.err <- 0.7
      })
      model({
        ka <- exp(tka + eta.ka)
        cl <- exp(tcl)
        v <- exp(tv)
        linCmt() ~ add(add.err)
      })
    }

    single <- nlmixr(one.cmt, theo_sd, "focei", list(print=0))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ the following etas are removed from the model since their initial estimates are zero: eta.ka
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> rxode2 3.0.2.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> calculating covariance matrix
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00 
#> done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress parHistData in nlmixr2 object, save 3960
#> ℹ initial model updated with final estimates, some zero etas are excluded from output
    print(single)
#> ── nlmixr² Population Only (outer: nlminb) ──
#> 
#>         OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> Pop 153.8231 404.4229 415.9541      -198.2115        19.97851        8.523242
#> 
#> ── Time (sec $time): ──
#> 
#>            setup optimize covariance table compress    other
#> elapsed 0.016396 0.001965   0.001966 0.023    0.009 1.378673
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>          Est.     SE %RSE Back-transformed(95%CI) BSV(SD) Shrink(SD)%
#> tka      0.44  0.213 48.3       1.55 (1.02, 2.35)                    
#> tcl     0.962  0.139 14.4       2.62 (1.99, 3.43)                    
#> tv       3.49 0.0881 2.52         32.9 (27.6, 39)                    
#> add.err  1.39                                1.39                    
#>  
#>   Covariance Type ($covMethod): r
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Information about run found ($runInfo):
#>    • gradient problems with initial estimate and covariance; see $scaleInfo 
#>   Censoring ($censInformation): No censoring
#>   Minimization message ($message):  
#>     relative convergence (4) 
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 11
#>   ID     TIME    DV IPRED   IRES  IWRES    ka    cl     v   tad dosenum
#>   <fct> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
#> 1 1      0     0.74  0     0.74   0.531  1.55  2.62  32.9  0          1
#> 2 1      0.25  2.84  3.10 -0.260 -0.186  1.55  2.62  32.9  0.25       1
#> 3 1      0.57  6.57  5.57  0.997  0.715  1.55  2.62  32.9  0.57       1
#> # ℹ 129 more rows

    single <- nlmixr(one.cmt, theo_sd, "bobyqa", list(print=0))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ the following etas are removed from the model since their initial estimates are zero: eta.ka
#> → pruning branches (`if`/`else`) of population log-likelihood model...
#> ✔ done
#> → loading llik model into symengine environment...
#> → finding duplicate expressions in population log-likelihood model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> |-----+---------------+-----------+-----------+-----------+-----------|
#> → calculating covariance
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> ✔ done
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → loading llik model into symengine environment...
#> → pruning branches (`if`/`else`) of llik full model...
#> ✔ done
#> → finding duplicate expressions in Llik EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in Llik EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling Llik EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress parHistData in nlmixr2 object, save 22064
#> ℹ initial model updated with final estimates, some zero etas are excluded from output
    print(single)
#> ── nlmixr² bobyqa ──
#> 
#>          OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> lPop 216.1508 466.7505 478.2818      -229.3753        19.98372        8.522614
#> 
#> ── Time (sec $time): ──
#> 
#>            setup table compress    other
#> elapsed 0.001546 0.022    0.009 2.480454
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>           Est.      SE  %RSE Back-transformed(95%CI) BSV(SD) Shrink(SD)%
#> tka     0.4403  0.2125 48.26    1.553 (1.024, 2.355)                    
#> tcl     0.9615  0.1387 14.43    2.616 (1.993, 3.433)                    
#> tv       3.492 0.08807 2.522    32.86 (27.65, 39.05)                    
#> add.err  1.375                                 1.375                    
#>  
#>   Covariance Type ($covMethod): r
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Censoring ($censInformation): No censoring
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 11
#>   ID     TIME    DV IPRED   IRES  IWRES    ka    cl     v   tad dosenum
#>   <fct> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
#> 1 1      0     0.74  0     0.74   0.538  1.55  2.62  32.9  0          1
#> 2 1      0.25  2.84  3.10 -0.260 -0.189  1.55  2.62  32.9  0.25       1
#> 3 1      0.57  6.57  5.57  0.996  0.724  1.55  2.62  32.9  0.57       1
#> # ℹ 129 more rows

    single <- nlmixr(one.cmt, theo_sd, "nls", list(print=0))
#> ℹ parameter labels from comments are typically ignored in non-interactive mode
#> ℹ Need to run with the source intact to parse comments
#> ℹ the following etas are removed from the model since their initial estimates are zero: eta.ka
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of nls model...
#> ✔ done
#> → calculate jacobian
#> → calculate ∂(f)/∂(θ)
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in nls gradient...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in nls gradient...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → finding duplicate expressions in nls pred-only...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> |-----+---------------+-----------+-----------+-----------|
#> → calculating covariance
#> ✔ done
#> → loading into symengine environment...
#> → pruning branches (`if`/`else`) of full model...
#> ✔ done
#> → finding duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → optimizing duplicate expressions in EBE model...
#> [====|====|====|====|====|====|====|====|====|====] 0:00:00
#> → compiling EBE model...
#> using C compiler: ‘gcc (Ubuntu 11.4.0-1ubuntu1~22.04) 11.4.0’
#> ✔ done
#> → Calculating residuals/tables
#> ✔ done
#> → compress origData in nlmixr2 object, save 5952
#> → compress parHistData in nlmixr2 object, save 2328
#> ℹ initial model updated with final estimates, some zero etas are excluded from output
    print(single)
#> ── nlmixr² nls with LM algorithm ──
#> 
#>         OBJF      AIC      BIC Log-likelihood Condition#(Cov) Condition#(Cor)
#> Pop 249.7119 500.3117 511.8429      -246.1558        21.78507        9.136566
#> 
#> ── Time (sec $time): ──
#> 
#>            setup table compress    other
#> elapsed 0.001683 0.034    0.008 2.299317
#> 
#> ── Population Parameters ($parFixed or $parFixedDf): ──
#> 
#>           Est.      SE  %RSE Back-transformed(95%CI) BSV(SD) Shrink(SD)%
#> tka     0.4403  0.1616  36.7    1.553 (1.132, 2.132)                    
#> tcl     0.9615   0.102 10.61    2.616 (2.141, 3.195)                    
#> tv       3.492 0.06569 1.881    32.86 (28.89, 37.38)                    
#> add.err  1.381                                 1.381                    
#>  
#>   Covariance Type ($covMethod): r (LM)
#>   No correlations in between subject variability (BSV) matrix
#>   Full BSV covariance ($omega) or correlation ($omegaR; diagonals=SDs) 
#>   Distribution stats (mean/skewness/kurtosis/p-value) available in $shrink 
#>   Censoring ($censInformation): No censoring
#>   Minimization message ($message):  
#>     Relative error in the sum of squares is at most `ftol'. 
#> 
#> ── Fit Data (object is a modified tibble): ──
#> # A tibble: 132 × 11
#>   ID     TIME    DV IPRED   IRES  IWRES    ka    cl     v   tad dosenum
#>   <fct> <dbl> <dbl> <dbl>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl>
#> 1 1      0     0.74  0     0.74   0.536  1.55  2.62  32.9  0          1
#> 2 1      0.25  2.84  3.10 -0.260 -0.189  1.55  2.62  32.9  0.25       1
#> 3 1      0.57  6.57  5.57  0.996  0.722  1.55  2.62  32.9  0.57       1
#> # ℹ 129 more rows

^{Created on 2024-11-11 with reprex v2.1.1}

[lirui0321]

Hi Matt,

There are a bunch of estimation methods listed for nlmixr2. Some are able to handle mixed effect model, some are not, if I understand them correctly. Is there a guidance about when to use which?

And one more question, should I do MU reference regardless of the estimation method? (NONMEM does not allow me to MU reference, if I only have a single subject)

[mattfidler]

BTW, I note that to update nlmixr2, I have to update R and Rtools from 4.3 to 4.4 as well, otherwise, new nlmixr2est would not load in old R.

This may be a windows issue, in our environment an older version of R runs. If you have issues (and post them), I could perhaps investigate.

[mattfidler]

There are a bunch of estimation methods listed for nlmixr2. Some are able to handle mixed effect model, some are not, if I understand them correctly. Is there a guidance about when to use which?

The mixed effect model estimation routines in nlmixr2 are:

• focei (and related methods)
• foce (= focei when residuals do not depend on the prediction as in the case of additive errors or lognormal errors)
• fo Poor man's mixed effect model (not the best, but could be used to quickly try out initial estimates and get better estimates...?). Can lead to different results when chaining results
• saem -- No gradients ever calculated, robust (but more likely to support over-paramterized models).
• nlme -- mixed effects models using R's built-in nlme toolbox

focei and family seem to perform better with low dimensional between subject variability mainly because it expands the ODE system by the number of compartments and number of between subject variability (Could change if I ever implement automatic-differentiation focei or adjoint sensitivity focei; classic focei in NONMEM and nlmixr2 uses forward sensitivities).

saem on the other hand works better with many between subject variabilities and over-parameterized models since it can search the space a bit better.

The rest of the methods are usually for population only estimation (though there are exceptions like simulation with rxSolve)

[mattfidler]

And one more question, should I do MU reference regardless of the estimation method? (NONMEM does not allow me to MU reference, if I only have a single subject)

Mu referencing regardless of estimation method is supported in nlmixr2, but what happens in mu referencing is dependent on the estimation method.

• saem combines the mu referenced expression to a single parameter and the optimizes for each individual. These values are used in a linear model to determine the population value, and the covariate values by a linear model. This is more stable. Note that 1 subject makes no sense in this context (and is why it isn't allowed in NONMEM).
• nlme uses the same method (and the linear models is likely what determines population values too).
• focei methods use the mu-referencing to determine theta drift and adjust the values, though in the future it may allow this linear model approach by a different estimation method.
• All other methods do not use mu-referencing for anything; There may be some stability enhancements if they are on a log-scale, though this is simply an effect of the parameterization instead of the linear models.

Hopefully this makes sense.

[lirui0321]

Thank you so much, Matt! Yes, that makes a lot of sense. NONMEM has an issue that MU referenced parameter must come with an ETA. It is not convenient if I want to switch from a population fitting to a single subject fitting, because I have to remove all MUs.

Rui