Error when fitting negative binomial model

billdenney commented 1 year ago

In the model below, I'm getting an error about the precision parameter, but I don't think I'm setting anything to -2147483648.

library(nlmixr2)
#> Loading required package: nlmixr2data

model_nb_linear_nlmixr <- function() {
  ini({
    ln0 <- log(3)
    sln <- -0.13
    lnbsize <- log(1.1)
  })
  model({
    nest <- exp(ln0 + sln*conc)
    nbsize <- exp(lnbsize)
    voc ~ dnbinomMu(nbsize, nest)
  })
}

d_model <-
  data.frame(
    ID = 0:9,
    TIME = 0:9,
    DV = rep(0:4, each = 2),
    conc = rep(0:4, each = 2)*exp(rnorm(10, mean = 0, sd = 0.1))
  )

nlmixr(object = model_nb_linear_nlmixr, data = d_model, est = "focei")
#> rxode2 2.0.13.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
#> Key: U: Unscaled Parameters; X: Back-transformed parameters; G: Gill difference gradient approximation
#> F: Forward difference gradient approximation
#> C: Central difference gradient approximation
#> M: Mixed forward and central difference gradient approximation
#> Unscaled parameters for Omegas=chol(solve(omega));
#> Diagonals are transformed, as specified by foceiControl(diagXform=)
#> |-----+---------------+-----------+-----------+-----------|
#> |    #| Objective Fun |       ln0 |       sln |   lnbsize |
#> |-----+---------------+-----------+-----------+-----------|
#> |    1|     198.74686 |     1.000 |    -1.000 |   -0.6332 |
#> |    U|     198.74686 |     1.099 |   -0.1300 |   0.09531 |
#> |    X|     198.74686 |     3.000 |   -0.1300 |     1.100 |
#> Error in .foceiFitInternal(.env) : 
#>   neg_binomial_2_lpmf: Precision parameter is -2147483648, but must be positive finite!
#> Restart 1
#> Error : focei$rxInv needs to be of class'rxSymInvCholEnv'
#> Error: focei$rxInv needs to be of class'rxSymInvCholEnv'

^{Created on 2023-08-29 with reprex v2.0.2}

mattfidler commented 1 year ago

This comes from stan and says you need to be more careful about your parameter estimates. I had a similar issue with a binomial distribution. I believe you will have to fix this in your model.

(ie, this is a user error and not a nlmixr2 error per se).

mattfidler commented 1 year ago

I believe the error is from here:

https://github.com/stan-dev/math/blob/bdf281f4e7f8034f47974d14dea7f09e600ac02a/stan/math/opencl/prim/neg_binomial_2_log_glm_lpmf.hpp#L159

mattfidler commented 1 year ago

I believe that NA equals -2147483648 though I cannot remember.

mattfidler commented 1 year ago

Indeed this comes from a NA integer converted to a double value

You can see the NA definition here:

https://github.com/wch/r-source/blob/2336bb0278ef721b9e22dc12d20770e85c6e48d1/src/main/arithmetic.c#L143

Which means –2147483648 is NA for an input integer.

See:

https://www.geeksforgeeks.org/int_max-int_min-cc-applications/

mattfidler commented 1 year ago

It could mean that integers that equal that value should be converted to Real NA values as well, though that is a large request

billdenney commented 1 year ago

Where would an NA come from in this model? All the values that I see should be positive since they are exponentiated.

mattfidler commented 1 year ago

There are embedded NAs in etTrans but nothing from an integer:

library(rxode2)
#> rxode2 2.0.13.9000 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

model_nb_linear_nlmixr <- function() {
  ini({
    ln0 <- log(3)
    sln <- -0.13
    lnbsize <- log(1.1)
  })
  model({
    nest <- exp(ln0 + sln*conc)
    nbsize <- exp(lnbsize)
    voc ~ dnbinomMu(nbsize, nest)
  })
}

d_model <-
  data.frame(
    ID = 0:9,
    TIME = 0:9,
    DV = rep(0:4, each = 2),
    conc = rep(0:4, each = 2)*exp(rnorm(10, mean = 0, sd = 0.1))
  )

etTrans(d_model, model_nb_linear_nlmixr)
#>    ID TIME EVID AMT II DV
#> 1   0    0    0  NA  0  0
#> 2   1    0    9  NA  0 NA
#> 3   1    1    0  NA  0  0
#> 4   2    0    9  NA  0 NA
#> 5   2    2    0  NA  0  1
#> 6   3    0    9  NA  0 NA
#> 7   3    3    0  NA  0  1
#> 8   4    0    9  NA  0 NA
#> 9   4    4    0  NA  0  2
#> 10  5    0    9  NA  0 NA
#> 11  5    5    0  NA  0  2
#> 12  6    0    9  NA  0 NA
#> 13  6    6    0  NA  0  3
#> 14  7    0    9  NA  0 NA
#> 15  7    7    0  NA  0  3
#> 16  8    0    9  NA  0 NA
#> 17  8    8    0  NA  0  4
#> 18  9    0    9  NA  0 NA
#> 19  9    9    0  NA  0  4
#> 
#> Covariates (non time-varying):
#>    ID      conc
#> 1   1 0.0000000
#> 2   2 0.0000000
#> 3   3 0.9393731
#> 4   4 1.2689009
#> 5   5 1.8698106
#> 6   6 1.6156174
#> 7   7 3.0705462
#> 8   8 2.6240829
#> 9   9 3.7391691
#> 10 10 3.5589855
#> 
#> Compartment translation:
#> [1] Compartment Number
#> <0 rows> (or 0-length row.names)

str(etTrans(d_model, model_nb_linear_nlmixr))
#> Classes 'rxEtTran' and 'data.frame': 19 obs. of  6 variables:
#>  $ ID  : Factor w/ 10 levels "0","1","2","3",..: 1 2 2 3 3 4 4 5 5 6 ...
#>  $ TIME: num  0 0 1 0 2 0 3 0 4 0 ...
#>  $ EVID: int  0 9 0 9 0 9 0 9 0 9 ...
#>  $ AMT : num  NA NA NA NA NA NA NA NA NA NA ...
#>  $ II  : num  0 0 0 0 0 0 0 0 0 0 ...
#>  $ DV  : num  0 NA 0 NA 1 NA 1 NA 2 NA ...

^{Created on 2023-08-29 with reprex v2.0.2}

billdenney commented 1 year ago

What is EVID = 9 (it doesn't show up here: https://nlmixr2.github.io/rxode2/articles/rxode2-event-types.html) and why is it inserted into the et table? Should it also have MDV = 1 associated with it so that the calculation is not performed?

mattfidler commented 1 year ago

EVID =9 does show up in the underlying true EVID structure here:

https://nlmixr2.github.io/rxode2/articles/rxode2-events-classic.html

mattfidler commented 1 year ago

MDV isn't carried along in the underlying translation.

mattfidler commented 1 year ago

Internally an EVID=9 is a non-observation event and makes sure the system is initialized to zero; EVID=9 should not be manually set. E

It shouldn't be affecting your problem. It has to come from an integer

mattfidler commented 1 year ago

If you simulate from the model, everything is NA because the parameters are not working for the negative binomial distribution as written.

library(rxode2)
#> rxode2 2.0.13.9000 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`
library(nlmixr2)
#> Loading required package: nlmixr2data

model_nb_linear_nlmixr <- function() {
  ini({
    ln0 <- log(3)
    sln <- -0.13
    lnbsize <- log(1.1)
  })
  model({
    nest <- exp(ln0 + sln*conc)
    nbsize <- exp(lnbsize)
    voc ~ dnbinomMu(nbsize, nest)
  })
}

d <-et(0:9) %>% et(id=1:10) %>% as.data.frame()
d$conc <- rep(0:4, each = 2)

ds <- rxSolve(model_nb_linear_nlmixr, d)

print(ds)
#> ── Solved rxode2 object ──
#> ── Parameters ($params): ──
#>         ln0         sln     lnbsize 
#>  1.09861229 -0.13000000  0.09531018 
#> ── Initial Conditions ($inits): ──
#> named numeric(0)
#> ── First part of data (object): ──
#> # A tibble: 100 × 7
#>      id  time  nest nbsize  pred  conc    DV
#>   <int> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
#> 1     1     0  3       1.1     8     0    NA
#> 2     1     1  3       1.1    11     0    NA
#> 3     1     2  2.63    1.1     0     1    NA
#> 4     1     3  2.63    1.1     1     1    NA
#> 5     1     4  2.31    1.1     4     2    NA
#> 6     1     5  2.31    1.1     2     2    NA
#> # ℹ 94 more rows

print(summary(ds))
#> ── Summary of Solved rxode2 object ──
#> ── Model ($model): ──
#> rxode2({
#>     param(ln0, sln, lnbsize, conc)
#>     nest = exp(ln0 + sln * conc)
#>     nbsize = exp(lnbsize)
#>     rx_yj_ ~ 182
#>     rx_lambda_ ~ 1
#>     rx_low_ ~ 0
#>     rx_hi_ ~ 1
#>     rx_r_ ~ 0
#>     rx_pred_ ~ llikNbinomMu(DV, nbsize, nest)
#>     sim = rxnbinomMu(nbsize, nest)
#>     cmt(voc)
#>     dvid(1)
#> }) 
#> ── Parameters (ds$params): ──
#>         ln0         sln     lnbsize 
#>  1.09861229 -0.13000000  0.09531018 
#> ── Initial Conditions (ds$inits): ──
#> named numeric(0)
#> ── Summary of data (object): ──
#>        id            time          nest           nbsize         pred      
#>  Min.   : 1.0   Min.   :0.0   Min.   :1.784   Min.   :1.1   Min.   : 0.00  
#>  1st Qu.: 3.0   1st Qu.:2.0   1st Qu.:2.031   1st Qu.:1.1   1st Qu.: 0.00  
#>  Median : 5.5   Median :4.5   Median :2.313   Median :1.1   Median : 1.00  
#>  Mean   : 5.5   Mean   :4.5   Mean   :2.352   Mean   :1.1   Mean   : 2.32  
#>  3rd Qu.: 8.0   3rd Qu.:7.0   3rd Qu.:2.634   3rd Qu.:1.1   3rd Qu.: 3.25  
#>  Max.   :10.0   Max.   :9.0   Max.   :3.000   Max.   :1.1   Max.   :14.00  
#>                                                                            
#>       conc         DV     
#>  Min.   :0   Min.   : NA  
#>  1st Qu.:1   1st Qu.: NA  
#>  Median :2   Median : NA  
#>  Mean   :2   Mean   :NaN  
#>  3rd Qu.:3   3rd Qu.: NA  
#>  Max.   :4   Max.   : NA  
#>              NA's   :100  
#>        id            time          nest           nbsize         pred      
#>  Min.   : 1.0   Min.   :0.0   Min.   :1.784   Min.   :1.1   Min.   : 0.00  
#>  1st Qu.: 3.0   1st Qu.:2.0   1st Qu.:2.031   1st Qu.:1.1   1st Qu.: 0.00  
#>  Median : 5.5   Median :4.5   Median :2.313   Median :1.1   Median : 1.00  
#>  Mean   : 5.5   Mean   :4.5   Mean   :2.352   Mean   :1.1   Mean   : 2.32  
#>  3rd Qu.: 8.0   3rd Qu.:7.0   3rd Qu.:2.634   3rd Qu.:1.1   3rd Qu.: 3.25  
#>  Max.   :10.0   Max.   :9.0   Max.   :3.000   Max.   :1.1   Max.   :14.00  
#>                                                                            
#>       conc         DV     
#>  Min.   :0   Min.   : NA  
#>  1st Qu.:1   1st Qu.: NA  
#>  Median :2   Median : NA  
#>  Mean   :2   Mean   :NaN  
#>  3rd Qu.:3   3rd Qu.: NA  
#>  Max.   :4   Max.   : NA  
#>              NA's   :100

^{Created on 2023-08-29 with reprex v2.0.2}

You could see if it not matching what you expect by looking at the source

https://github.com/nlmixr2/rxode2ll/blob/main/src/llikNbinom2.cpp

And the description of the distribution from stan:

https://mc-stan.org/docs/2_20/functions-reference/nbalt.html

It could be that I'm implementing this incorrectly, but I'm only pulling things from stan right now for likelihoods.

billdenney commented 1 year ago

When I tested my parameters in R with what I think is the same parameterization, I got what I expected out:

> dnbinom(0:10, size = 1.1, mu = 3)
 [1] 0.23521754 0.18932144 0.14545428 0.10997762 0.08248322 0.06156065 0.04579511 0.03398731 0.02517963 0.01862883 0.01376716

mattfidler commented 1 year ago

dnbinom is a different distribution since it is from R. You need to test the stan parameterization, or the underlying function:

library(rxode2)
#> rxode2 2.0.13.9000 using 4 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

model_nb_linear_nlmixr <- function() {
  ini({
    ln0 <- log(3)
    sln <- -0.13
    lnbsize <- log(1.1)
  })
  model({
    nest <- exp(ln0 + sln*conc)
    nbsize <- exp(lnbsize)
    voc ~ dnbinomMu(nbsize, nest)
  })
}

f <- model_nb_linear_nlmixr()

summary(f$simulationModel)
#> rxode2 2.0.13.9000 model named rx_a2d71a872d72c8620009fbf2ee6c6790_x6 model (✔ ready). 
#> DLL: C:\Users\fidlema3\AppData\Local\Temp\1\RTMPCA~1\rxode2\RX_A2D~1.RXD/rx_a2d71a872d72c8620009fbf2ee6c6790_x64.dll
#> NULL
#> 
#> Calculated Variables:
#> [1] "nest"   "nbsize" "sim"   
#> ── rxode2 Model Syntax ──
#> rxode2({
#>     param(ln0, sln, lnbsize, conc)
#>     nest = exp(ln0 + sln * conc)
#>     nbsize = exp(lnbsize)
#>     rx_yj_ ~ 182
#>     rx_lambda_ ~ 1
#>     rx_low_ ~ 0
#>     rx_hi_ ~ 1
#>     rx_r_ ~ 0
#>     rx_pred_ ~ llikNbinomMu(DV, nbsize, nest)
#>     sim = rxnbinomMu(nbsize, nest)
#>     cmt(voc)
#>     dvid(1)
#> })

^{Created on 2023-08-30 with reprex v2.0.2}

mattfidler commented 1 year ago

It is likely because size is not an integer, which seems to be required.

rxode2random::rxnbinomMu(1.1, 3, n=100)
#> Error in rxode2random::rxnbinomMu(1.1, 3, n = 100): Assertion on 'size' failed: Must be of type 'count', not 'double'.

rxode2random::rxnbinomMu(1, 3, n=100)
#>   [1] 0 1 1 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 1 1 0 0 0 0 0 0 0 1 0 0
#>  [38] 0 1 0 1 0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 1 0 0 1 1 0 0 1 0 0 0 0 1 0 1 0 0 0
#>  [75] 1 0 0 0 0 0 1 0 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 1

^{Created on 2023-08-30 with reprex v2.0.2}

nlmixr2 / rxode2ll

Error when fitting negative binomial model #6