gnlmm(
  llik,
  data,
  inits,
  syspar = NULL,
  system = NULL,
  diag.xform = c("sqrt", "log", "identity"),
  ...,
  control = list()
)

gnlmm2(
  llik,
  data,
  inits,
  syspar = NULL,
  system = NULL,
  diag.xform = c("sqrt", "log", "identity"),
  ...,
  control = list()
)

if (FALSE) {
llik <- function() {
  lp <- THETA[1] * x1 + THETA[2] * x2 + (x1 + x2 * THETA[3]) * ETA[1]
  p <- pnorm(lp)
  dbinom(x, m, p, log = TRUE)
}
inits <- list(THTA = c(1, 1, 1), OMGA = list(ETA[1] ~ 1))

try(gnlmm(llik, rats, inits, control = list(nAQD = 1)))

llik <- function() {
  if (group == 1) {
    lp <- THETA[1] + THETA[2] * logtstd + ETA[1]
  } else {
    lp <- THETA[3] + THETA[4] * logtstd + ETA[1]
  }
  lam <- exp(lp)
  dpois(y, lam, log = TRUE)
}
inits <- list(THTA = c(1, 1, 1, 1), OMGA = list(ETA[1] ~ 1))

fit <- try(gnlmm(llik, pump, inits,
  control = list(
    reltol.outer = 1e-4,
    optim.outer = "nmsimplex",
    nAQD = 5
  )
))



ode <- "
d/dt(depot) =-KA*depot;
d/dt(centr) = KA*depot - KE*centr;
"
sys1 <- RxODE(ode)

pars <- function() {
  CL <- exp(THETA[1] + ETA[1]) # ; if (CL>100) CL=100
  KA <- exp(THETA[2] + ETA[2]) # ; if (KA>20) KA=20
  KE <- exp(THETA[3])
  V <- CL / KE
  sig2 <- exp(THETA[4])
}
llik <- function() {
  pred <- centr / V
  dnorm(DV, pred, sd = sqrt(sig2), log = TRUE)
}
inits <- list(THTA = c(-3.22, 0.47, -2.45, 0))
inits$OMGA <- list(ETA[1]+ETA[2]~c(.027, .01, .37))

theo <- theo_md

fit <- try(gnlmm(llik, theo, inits, pars, sys1,
  control = list(trace = TRUE, nAQD = 1)
))

fit2 <- try(gnlmm2(llik, theo, inits, pars, sys1,
  control = list(trace = TRUE, nAQD = 1)
))

if (inherits(fit, "gnlmm.fit")) {
 cv <- calcCov(fit)
 cbind(fit$par[fit$nsplt == 1], sqrt(diag(cv)))
}
}