'nlme' support

[mattfidler]

While linear mixed effect models are supported, nonlinear mixed effects models like nlme are not supported at this time. I'm unsure this is in scope or out of scope for this project.

[strengejacke]

Related #66
It's definitely inside the project's scope.

[mattfidler]

I would add nlmixr support too. But I'm unsure if all the stats used are always used by nonlinear mixed effects models.

What are the generics you are trying to extract? There doesn't seem to be a vignette like broom or broom.mixed has.

[strengejacke]

Most important functions are

• find_formula()
• get_paramters() and find_parameters()
• n_obs()
• get_data()

[mattfidler]

What is find_formula used for? In nlmixr there is not a formula, but a model.

[strengejacke]

find_formula() is a generic replacement of stats:: formula(), but splitting the parts (if any) into different list elements (i.e. formula for fixed effects, for random effects, for zero inflation...). It's the base for many other functions like find_terms() etc.

[DominiqueMakowski]

@mattfidler if you can post here a reproducible example with simple models + how to access the above information for them (formula, data etc), we can then easily implement its support :)

[mattfidler]

I am more familiar with nlmixr, give me a few weeks and I will give a simple example

[femiguez]

Not sure if this is useful, but in package nlraa I have a function 'var_cov' which extracts the variance-covariance matrix of any object produced by the nlme package. Notice that the function 'nlme::getVarCov' is limited and it fails for a variety of cases. Let me know if I can help with this.

[strengejacke]

Not sure if this is useful, but in package nlraa I have a function 'var_cov' which extracts the variance-covariance matrix of any object produced by the nlme package. Notice that the function 'nlme::getVarCov' is limited and it fails for a variety of cases. Let me know if I can help with this.

That sounds great, and it sounds like it could be useful here #206 (for insight::get_variance())

[femiguez]

I also tried ggeffects a while ago and noticed that it does not handle nlme objects (at least as I expected). The function 'nlraa::predict_nlme' produced confidence bands and could be called by 'ggeffects'

[mattfidler]

For nlmixr these components are n the code below:

library(nlmixr)

## simple nlmixr model
one.cmt <- function() {
  ini({
    ## You may label each parameter with a comment
    tka <- 0.45 # Log Ka
    tcl <- log(c(0, 2.7, 100)) # Log Cl
    ## This works with interactive models
    ## You may also label the preceding line with label("label text")
    tv <- 3.45; label("log V")
    ## the label("Label name") works with all models
    eta.ka ~ 0.6
    eta.cl ~ 0.3
    eta.v ~ 0.1
    add.sd <- 0.7
  })
  model({
    ka <- exp(tka + eta.ka)
    cl <- exp(tcl + eta.cl)
    v <- exp(tv + eta.v)
    linCmt() ~ add(add.sd)
  })
}

fit <- nlmixr(one.cmt, theo_sd, list(print=0), est="focei")
#> ℹ 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 
#> done
#> Calculating residuals/tables
#> done
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : initial ETAs were nudged; (can control by foceiControl(etaNudge=.,
#> etaNudge2=))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : ETAs were reset to zero during optimization; (Can control by
#> foceiControl(resetEtaP=.))
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod
#> = FALSE, : last objective function was not at minimum, possible problems in
#> optimization
#> Warning in (function (uif, data, est = NULL, control = list(), ..., sum.prod =
#> FALSE, : gradient problems with initial estimate and covariance; see $scaleInfo

library(broom.mixed)
#> Registered S3 method overwritten by 'broom.mixed':
#>   method      from 
#>   tidy.gamlss broom

# `get_parameters()`

# Not clear from the documentation exactly what this means:
# https://easystats.github.io/insight/reference/get_parameters.html

tidy(fit)
#> # A tibble: 7 x 7
#>   effect   group         term       estimate std.error statistic    p.value
#>   <chr>    <chr>         <chr>         <dbl>     <dbl>     <dbl>      <dbl>
#> 1 fixed    <NA>          tka           0.464    0.195       2.38  9.51e-  3
#> 2 fixed    <NA>          tcl           1.01     0.0751     13.5   2.04e- 26
#> 3 fixed    <NA>          tv            3.46     0.0437     79.3   5.01e-109
#> 4 ran_pars ID            sd__eta.ka    0.635   NA          NA    NA        
#> 5 ran_pars ID            sd__eta.cl    0.263   NA          NA    NA        
#> 6 ran_pars ID            sd__eta.v     0.138   NA          NA    NA        
#> 7 ran_pars Residual(add) add.sd        0.695   NA          NA    NA

fit$parFixed
#>        Parameter  Est.     SE %RSE Back-transformed(95%CI) BSV(CV%) Shrink(SD)%
#> tka       Log Ka 0.464  0.195 42.1       1.59 (1.08, 2.33)     70.4      1.81%<
#> tcl       Log Cl  1.01 0.0751 7.42       2.75 (2.37, 3.19)     26.8      3.87%<
#> tv         log V  3.46 0.0437 1.26       31.8 (29.2, 34.6)     13.9      10.3%<
#> add.sd           0.695                               0.695

fit$parFixedDf
#>         Estimate         SE      %RSE Back-transformed  CI Lower  CI Upper
#> tka    0.4637611 0.19520029 42.090698        1.5900431  1.084561  2.331115
#> tcl    1.0118846 0.07508740  7.420549        2.7507803  2.374332  3.186915
#> tv     3.4595102 0.04365043  1.261752       31.8013977 29.193818 34.641886
#> add.sd 0.6945831         NA        NA        0.6945831        NA        NA
#>        BSV(CV%) Shrink(SD)%
#> tka    70.42453    1.807129
#> tcl    26.75077    3.872229
#> tv     13.89715   10.303655
#> add.sd       NA          NA


# `n_obs()`

nobs(fit)
#> [1] 132


# `get_data()`

head(getData(fit))
#>   ID TIME    DV     AMT EVID CMT   WT
#> 1  1 0.00  0.00 319.992  101   1 79.6
#> 2  1 0.00  0.74   0.000    0   2 79.6
#> 3  1 0.25  2.84   0.000    0   2 79.6
#> 4  1 0.57  6.57   0.000    0   2 79.6
#> 5  1 1.12 10.50   0.000    0   2 79.6
#> 6  1 2.02  9.66   0.000    0   2 79.6


# find_formula()
# Closest  is
fit$uif
#> ▂▂ RxODE-based 1-compartment model with first-order absorption ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂ 
#> ── Initialization: ───────────────────────────────────────────────────────────── 
#> Fixed Effects ($theta): 
#>       tka       tcl        tv 
#> 0.4500000 0.9932518 3.4500000 
#> 
#> Omega ($omega): 
#>        eta.ka eta.cl eta.v
#> eta.ka    0.6    0.0   0.0
#> eta.cl    0.0    0.3   0.0
#> eta.v     0.0    0.0   0.1
#> ── μ-referencing ($muRefTable): ──────────────────────────────────────────────── 
#>                              ┌─────────┬─────────┐
#>                              │ theta   │ eta     │
#>                              ├─────────┼─────────┤
#>                              │ tka     │ eta.ka  │
#>                              ├─────────┼─────────┤
#>                              │ tcl     │ eta.cl  │
#>                              ├─────────┼─────────┤
#>                              │ tv      │ eta.v   │
#>                              └─────────┴─────────┘
#> 
#> ── Model: ────────────────────────────────────────────────────────────────────── 
#>     ka <- exp(tka + eta.ka)
#>     cl <- exp(tcl + eta.cl)
#>     v <- exp(tv + eta.v)
#>     linCmt() ~ add(add.sd) 
#> ▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂▂

^{Created on 2021-04-14 by the reprex package (v2.0.0)}

[mattfidler]

If you choose to use broom variants, there is more documentation here:

https://nlmixrdevelopment.github.io/nlmixr/articles/broom.html

[emphillips18]

Are there any updates on getting nlme support? (I apologize if this is not the space to ask, I am new to github =P)

[femiguez]

@emphillips18 Is there anything you need not available in package 'nlraa'?

[emphillips18]

I've never tried it. I'll take a look at it and get back to you. Thanks!

…

On Thu, Aug 17, 2023, 04:02 Fernando Miguez ***@***.***> wrote: @emphillips18 <https://github.com/emphillips18> Is there anything you need not available in package 'nlraa'? — Reply to this email directly, view it on GitHub <#309 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AMPA2HOQ4JRA7QS7PBWIINLXVX23LANCNFSM4YKC5VDQ> . You are receiving this because you were mentioned.Message ID: ***@***.***>