troubleshooting tips - Zero `FIM` determinant about poped HOT 13 CLOSED

Hi Vijay,

This generally happens when one parameter is not identifiable in the model. You can see where the problem is by looking at the FIM and finding the row/column that has all zeros in it. The order of the parameters in the FIM is c(pop,d,sigma). If you provide the ff fg and ferror functions then I can give more concrete advice.

Best regards,
Andy

On 30 May 2016, at 7:45 , Vijay Ivaturi [email protected] wrote:

hi Andy,
I am not sure how to diagnose the condition where the evaluate.fim returns a 0 determinant. I know that my ff_file returns a valid y vector (checked by running the setup outside a function setting). The poped.db seems to be have been set up correctly as far as the design goes, this I can confirm by plotting the predictions using the plot_model_prediction which returns a population prediction at the specified design times. But when I run an evaluate.fim on this poped.db the return value is a matrix of zero values and hence the determinant of zero.

I want to go about diagnosing why this happens, and was wondering if you could provide some tips. Here is my create.poped.db function. Just a note that A and B below are parameters that KA is derived from.

poped.db <- create.poped.database(ff_file="ff",
fg_file="sfg",
fError_file="feps",
bpop=c(CL=10,VC=60.1,Q=2,VP=30,TVLAG=0.2,A=0.256,B=0.284),
notfixed_bpop=c(1,1,1,1,0,0,0),
d=c(CL=0.08,VC=0.1),
sigma=c(0.06,6.48),
notfixed_sigma=c(0,0),
m=1,
groupsize=6,
xt=c(3,5,8,10),
minxt=c(3,4,7,9),
maxxt=c(4,7,9,10),
bUseGrouped_xt=0,
a = c(40000,70),
maxa=c(40000,70),
mina=c(40000,70))
—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub #2, or mute the thread https://github.com/notifications/unsubscribe/AFLmZP83gHZhsKwXbpJ3xoBv2-fzBHPwks5qGnlugaJpZM4IpgX_.

from poped.

Hmm.. I am getting zeros for all rows and columns. Let me re-check and get back to you, probably in the evening as I am out now :/

from poped.

Here are the ff, fg and feps functions - can't seem to figure it out. I get all columns as zero in my FIM

ff

ff <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)  
  dose_times <- 0
  time <- sort(unique(c(times_xt,dose_times)))
  is.dose <- time %in% dose_times

  data <- 
    dplyr::data_frame(ID = 1,
                      time = time,
                      amt = ifelse(is.dose,p[["DOSE"]], 0), 
                      cmt = ifelse(is.dose, 1, 0), 
                      evid = cmt,
                      CL = p[["CL"]], VC = p[["VC"]], Q = p[["Q"]],
                      VP = p[["VP"]], TVLAG = p[["TVLAG"]], A = p[["A"]], B = p[["B"]])

  out <- mrgsim(mod, data=data,param=as.list(p),recsort=4)

  y <-  out$CP

  y <- y[match(times_xt,out$time)]

  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
}

fg and feps

## -- parameter definition function 
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                VC=bpop[2]*exp(b[2]),
                Q=bpop[3],
                VP=bpop[4],
                TVLAG=bpop[5],
                A=bpop[6],
                B=bpop[7],
                DOSE=a[1],
                WT=a[2])
  return( parameters ) 
}
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional  
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])+epsi[,2]
  return(list( y= y,poped.db =poped.db )) 
}

poped.db

poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=10,VC=60.1,Q=2,VP=30,TVLAG=0.2,A=0.256,B=0.284), 
                                  notfixed_bpop=c(1,1,1,1,0,0,0),
                                  d=c(CL=0.08,VC=0.1), 
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=1,
                                  groupsize=6,
                                  xt=c(3,5,8,10),
                                  minxt=c(3,4,7,9),
                                  maxxt=c(4,7,9,10),
                                  bUseGrouped_xt=0,
                                  a = c(40000,70),
                                  maxa=c(40000,70),
                                  mina=c(40000,70))

from poped.

still trying to diagnose and going through the functions one at a time. I see that for some reason the bpop matrix in my poped.db looks like this below. Based on your statement above that " The order of the parameters in the FIM is c(pop,d,sigma)", I wonder if the bpop is being read as zero? Why would I see it this way, when I clearly specified my parameters as bpop in the create.poped.db function?

UPDATE: never mind - realized that poped.db.parameters follows the c(pop,d,sigma) pattern and not the poped.db.parameters$bpop. search continues..

from poped.

I saw something like this when I was coding the model in mrgsolve as well. Trying to recall what fixed the issue.

I'll take a look this morning.

One note you shouldn't have to include parameters on the data set; if I recall, that parameter set is just a single set of values. Updating the base parameter list (mrgsim(...., param = as.list(p), ...) should be sufficient. But let me know if as.list(p) isn't simple like that.

Can you post the model code too?

from poped.

Thanks, Kyle - here is the model code.

$CMT DEPOT CENT PERI

$SET end=168, delta=0.1
$PARAM TVCL=20.9, TVVC=70.1,TVQ=1, TVVP=27.5, TVLAG=0.175, WT=70,
TVA = 0.256 , TVB = 0.284

$OMEGA
prefix="E", labels = s(CL,VC,A,B)
0 0 0 0 

$MAIN if(EVID==1) double TOD=TIME;
double TAD = 0;
if(TIME>0) TAD = TIME - TOD;

double ASCL = pow(WT/70,0.75);
double ASV = pow(WT/70,1);

double CL=TVCL*ASCL*exp(ECL);

double VC=TVVC*ASV*exp(EVC);

double Q=TVQ*ASCL;

double VP=TVVP*ASV;

double A = TVA*exp(EA);
double B = TVB*exp(EB);

double TADCO = 24;
if(TAD<24) TADCO=TAD;

double KA = A + B*TADCO;

_ALAG(1) = TVLAG;

_F(2) = 1;


$ODE
double CP =  (CENT/VC);
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = -(CL/VC)*CENT -(Q/VC)*CENT + (Q/VP)*PERI + KA*DEPOT;
dxdt_PERI = (Q/VC)*CENT - (Q/VP)*PERI;

$CAPTURE CP
'

from poped.

@vjd

First:
The parameters that get passed into ff have names CL, VC, etc ... but the parameters in the model code are TVCL, TVVC, etc. So most of the parameters in that p list you pass to update won't actually update: the names don't match. This is why you are getting zero FIM.

Second:
A related issue is where the covariate model should happen. You have it in the mrgsolve model, but it sort of seems to me like it should be happening in the parameter definition function? I don't know enough about PopED to say one way or another. In that case, it seems to me like the model needs to be written in terms of basic PK parameters (CL, VC, Q etc) and let PopED realize the covariate-adjusted individual parameters (in that sfg function).

If it was just the first issue, we could take this off the PopED issue tracker. But I'd like to get some resolution / clarification on the second issue.

If you have:

a = c(40000,70),
maxa=c(40000,70),
mina=c(40000,70))

in the data base, does WT even matter?

from poped.

This gets you non-zero FIM.

library(mrgsolve) #latest from github
library(PopED)
library(PKPDsim)
library(deSolve)
library(ggplot2)
library(dplyr)



code_mod <- '

$CMT DEPOT CENT PERI

$PARAM 
CL=20.9, VC=70.1,Q=1, VP=27.5, TVLAG=0.175, WT=70,
KA = 0.256

$MAIN 
ALAG_DEPOT = TVLAG;
F_CENT = 1;

$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = -(CL/VC)*CENT -(Q/VC)*CENT + (Q/VP)*PERI + KA*DEPOT;
dxdt_PERI = (Q/VC)*CENT - (Q/VP)*PERI;

$TABLE
table(CP) = CENT/VC;

'

mod <- mread(code=code_mod, model="ex2_mrg_poped")

data <- ev(amt=c(4,6,8,15), cmt=1) %>% as.data.frame %>% as.tbl %>% mutate(ID=1:4, dose = amt*10,amt=amt*10000)
out <- mod %>% data_set(data) %>% Req(CP) %>% carry.out(dose) %>% mrgsim(end=48, delta=0.1)

plot(out)

fast <- TRUE 

iNumSimulations <- ifelse(fast,5,100)
EAStepSize <- ifelse(fast,40,1)
rsit <- ifelse(fast,3,300)
sgit <- ifelse(fast,3,150)
ls_step_size <- ifelse(fast,3,50)
iter_max <- ifelse(fast,1,10)

ff <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)  
  dose_times <- 0
  time <- sort(unique(c(times_xt,dose_times)))
  is.dose <- time %in% dose_times

  data <- 
    dplyr::data_frame(ID = 1,
                      time = time,
                      amt = ifelse(is.dose,p[["DOSE"]], 0), 
                      cmt = ifelse(is.dose, 1, 0), 
                      evid = cmt)

  out <- mrgsim(mod,end=-1, data=data,param=as.list(p),recsort=4)

  y <-  out$CP

  y <- y[match(times_xt,out$time)]

  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
}

## Define other functions


## -- parameter definition function 
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                VC=bpop[2]*exp(b[2]),
                Q=bpop[3],
                VP=bpop[4],
                TVLAG=bpop[5],
                DOSE=a[1],
                WT=a[2])
  return( parameters ) 
}
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional  
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db)) 
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])+epsi[,2]
  return(list( y= y,poped.db =poped.db )) 
}


poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175), 
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1), 
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=1,
                                  groupsize=6,
                                  xt=c(3,5,8,18),
                                  minxt=c(3,4,7,9),
                                  maxxt=c(4,7,9,24),
                                  bUseGrouped_xt=0,
                                  a = c(40000,70),
                                  maxa=c(40000,70),
                                  mina=c(40000,70))


system.time(p1 <- plot_model_prediction(poped.db))
p1


system.time(p2 <- plot_model_prediction(poped.db,IPRED=TRUE,IPRED.lines=FALSE, IPRED.lines.pctls=FALSE,
                      DV=TRUE, DV.lines=FALSE, DV.points=FALSE, separate.groups=FALSE, PRED=TRUE,
                      sample.times = TRUE, sample.times.IPRED = FALSE, sample.times.DV = FALSE))
p2


### Evaluate FIM
FIM <- evaluate.fim(poped.db, fim.calc.type = 1, deriv.type = 1) 
det(FIM)
get_rse(FIM, poped.db)

But get_rse has very large values and I get this error when running the optimization:

*******************************************
Stopping criteria testing
*******************************************
Efficiency of design (start of loop / end of loop) =  NaN 
Efficiency stopping criteria (lower limit) =  1.001 
Error in if (eff <= eff_crit) stop_crit <- TRUE : 
  missing value where TRUE/FALSE needed

from poped.

Thanks @kylebmetrum. A couple of changes. You had KA defined as a parameter in the model, but if was not present in sfg or poped.db. I went ahead made these changes, but as you pointed out, I am also getting the really large values for get_rse. Don't get the message that you have.
Any insights, @andrewhooker ?

from poped.

@vjd Seems like you can get rse down by adding a couple of points farther out in time (it's only allowed to go out to 24 hours right now). Also, I wonder if you need to either get more intensive on the single dose or do multiple doses to get information on VP. But if @andrewhooker chimes in on this, listen to him, not me.

This was successful run for me.

poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175), 
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1), 
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=1,
                                  groupsize=20,
                                  xt=c(      1,3,5,8,28,96),
                                  minxt=c(0.01,3,4,7,24,72),
                                  maxxt=c(2,   4,7,9,36,120),
                                  bUseGrouped_xt=0,
                                  a = c(40000,70),
                                  maxa=c(40000,70),
                                  mina=c(40000,70))

FINAL RESULTS
Optimized Sampling Schedule
0.4347      3      7      9     24     72

OFV = 3.64853

Efficiency [typically: (OFV_final/OFV_initial)^(1/npar)]: 1.05654

Expected parameter 
relative standard error (%RSE):
   Parameter   Values    RSE_0      RSE
     bpop[1]    20.90     7.98     7.39
     bpop[2]    70.10     8.85     8.65
     bpop[3]     1.00    70.32    60.89
     bpop[4]    27.50   209.68   142.26
      D[1,1]     0.08    41.43    37.83
      D[2,2]     0.10    46.21    43.92

Total running time: 128.334 seconds
> get_rse(output$FIM,output$poped.db)
   bpop[1]    bpop[2]    bpop[3]    bpop[4]     D[1,1]     D[2,2] 
  7.394085   8.650731  60.892304 142.260853  37.833258  43.924750 

from poped.

Hi

I think the issue is that you have sparse sampling (4 samples per individual) and only one group (so all individuals have the sample sampling time), plus 4 + 2 parameters to estimate. You need either more samples per individual or groups of individuals with different samples, plus potentially a later sampling time frame. In general, if your initial design is so bad that you have thousand’s of percent RSE predicted, and you have a relatively restricted design (you had only a small range of allowed time points per sample) then you can’t expect to get good results from an optimization.

I attach code that optimizes 2 groups with 3 ids each and sample times between 0 and 120, and get:

FINAL RESULTS
Optimized Sampling Schedule
Group 1 : 0.5211  1.077   24.7  74.41
Group 2 :  1.667  15.91  31.99  76.54

OFV = -4.87805

Efficiency [typically: (OFV_final/OFV_initial)^(1/npar)]: NaN

D-Efficiency [(det(FIM_final)/det(FIM_initial))^(1/npar)]: 1.58166

Expected parameter
relative standard error (%RSE):
   Parameter   Values   RSE_0     RSE
     bpop[1]    20.90    15.0    16.8
     bpop[2]    70.10    19.1    16.0
     bpop[3]     1.00   191.3   107.9
     bpop[4]    27.50   376.3   260.6
      D[1,1]     0.08    84.2    82.2
      D[2,2]     0.10    97.3    85.1

Total running time: 208.122 seconds

Still poor estimates of pop[3 and 4] but better. Negative efficiency is just a consequence of having negative objective functions (very low information).

Code below

#devtools::install_github("metrumresearchgroup/mrgsolve", subdir="rdev")
#devtools::install_github("ronkeizer/PKPDsim")

packageVersion("PopED")
packageVersion("mrgsolve")
packageVersion("PKPDsim")

library(mrgsolve)
library(PopED)
library(PKPDsim)
library(deSolve)
library(ggplot2)
library(dplyr)



code_mod <- '

$CMT DEPOT CENT PERI

$PARAM
CL=20.9, VC=70.1,Q=1, VP=27.5, TVLAG=0.175, WT=70,
KA = 0.256

$MAIN
ALAG_DEPOT = TVLAG;
F_CENT = 1;

$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = -(CL/VC)*CENT -(Q/VC)*CENT + (Q/VP)*PERI + KA*DEPOT;
dxdt_PERI = (Q/VC)*CENT - (Q/VP)*PERI;

$TABLE
table(CP) = CENT/VC;

'

mod <- mread(code=code_mod, model="ex2_mrg_poped")

data <- ev(amt=c(4,6,8,15), cmt=1) %>% as.data.frame %>% as.tbl %>% mutate(ID=1:4, dose = amt*10,amt=amt*10000)
out <- mod %>% data_set(data) %>% Req(CP) %>% carry.out(dose) %>% mrgsim(end=48, delta=0.1)

plot(out)

ff <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)
  dose_times <- 0
  time <- sort(unique(c(times_xt,dose_times)))
  is.dose <- time %in% dose_times

  data <-
    dplyr::data_frame(ID = 1,
                      time = time,
                      amt = ifelse(is.dose,p[["DOSE"]], 0),
                      cmt = ifelse(is.dose, 1, 0),
                      evid = cmt)

  out <- mrgsim(mod,end=-1, data=data,param=as.list(p),recsort=4)

  y <-  out$CP

  y <- y[match(times_xt,out$time)]

  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
}

## Define other functions


## -- parameter definition function
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                VC=bpop[2]*exp(b[2]),
                Q=bpop[3],
                VP=bpop[4],
                TVLAG=bpop[5],
                DOSE=a[1],
                WT=a[2])
  return( parameters )
}
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])+epsi[,2]
  return(list( y= y,poped.db =poped.db ))
}


poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175),
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1),
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=1,
                                  groupsize=6,
                                  xt=c(3,5,8,18),
                                  minxt=c(3,4,7,9),
                                  maxxt=c(4,7,9,24),
                                  bUseGrouped_xt=0,
                                  a = c(40000,70),
                                  maxa=c(40000,70),
                                  mina=c(40000,70))


system.time(p1 <- plot_model_prediction(poped.db))
p1


system.time(p2 <- plot_model_prediction(poped.db,IPRED=TRUE,IPRED.lines=FALSE, IPRED.lines.pctls=FALSE,
                                        DV=TRUE, DV.lines=FALSE, DV.points=FALSE, separate.groups=FALSE, PRED=TRUE,
                                        sample.times = TRUE, sample.times.IPRED = FALSE, sample.times.DV = FALSE))
p2


### Evaluate FIM
FIM <- evaluate.fim(poped.db)
FIM
det(FIM)
get_rse(FIM, poped.db)

## here we see that lots of parameters are poorly estimated
## not surprising, 1 group, 4 samples, and 4 bpop and 2 d parameters to estimate
poped_db_2 <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175),
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1),
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=2,
                                  groupsize=3,
                                  xt=list(c(3,5,8,18),c(1,6,50,120)),
                                  minxt=0,
                                  maxxt=120,
                                  a = c(40000,70))
FIM <- evaluate.fim(poped_db_2)
FIM
det(FIM)
get_rse(FIM, poped_db_2)


# run only methods that can be parallelized to make things faster
poped_optim(poped_db_2, opt_xt = T,parallel = T, method=c("ARS"))  # approximate
#poped_optim(poped_db_2, opt_xt = T,parallel = T, method=c("ARS","GA","LS")) # more accurate

from poped.

the design space was indeed poor. In the spirit of trying out, I gave a few sample time points without much consideration before getting into the actual sample time optimization in different pediatric age groups that I was attempting (each group with a different median body weight). Your implementation is a good start for me to play around further. Thank you @kylebmetrum for helping with the mrgsolve aspect.

from poped.

Hi

I think the issue is that you have sparse sampling (4 samples per individual) and only one group (so all individuals have the sample sampling time), plus 4 + 2 parameters to estimate. You need either more samples per individual or groups of individuals with different samples, plus potentially a later sampling time frame. In general, if your initial design is so bad that you have thousand’s of percent RSE predicted, and you have a relatively restricted design (you had only a small range of allowed time points per sample) then you can’t expect to get good results from an optimization.

I attach code that optimizes 2 groups with 3 ids each and sample times between 0 and 120, and get:

FINAL RESULTS
Optimized Sampling Schedule
Group 1 : 0.5211  1.077   24.7  74.41
Group 2 :  1.667  15.91  31.99  76.54

OFV = -4.87805

Efficiency [typically: (OFV_final/OFV_initial)^(1/npar)]: NaN

D-Efficiency [(det(FIM_final)/det(FIM_initial))^(1/npar)]: 1.58166

Expected parameter
relative standard error (%RSE):
   Parameter   Values   RSE_0     RSE
     bpop[1]    20.90    15.0    16.8
     bpop[2]    70.10    19.1    16.0
     bpop[3]     1.00   191.3   107.9
     bpop[4]    27.50   376.3   260.6
      D[1,1]     0.08    84.2    82.2
      D[2,2]     0.10    97.3    85.1

Total running time: 208.122 seconds

Still poor estimates of pop[3 and 4] but better. Negative efficiency is just a consequence of having negative objective functions (very low information).

Code below

#devtools::install_github("metrumresearchgroup/mrgsolve", subdir="rdev")
#devtools::install_github("ronkeizer/PKPDsim")

packageVersion("PopED")
packageVersion("mrgsolve")
packageVersion("PKPDsim")

library(mrgsolve)
library(PopED)
library(PKPDsim)
library(deSolve)
library(ggplot2)
library(dplyr)



code_mod <- '

$CMT DEPOT CENT PERI

$PARAM
CL=20.9, VC=70.1,Q=1, VP=27.5, TVLAG=0.175, WT=70,
KA = 0.256

$MAIN
ALAG_DEPOT = TVLAG;
F_CENT = 1;

$ODE
dxdt_DEPOT = -KA*DEPOT;
dxdt_CENT = -(CL/VC)*CENT -(Q/VC)*CENT + (Q/VP)*PERI + KA*DEPOT;
dxdt_PERI = (Q/VC)*CENT - (Q/VP)*PERI;

$TABLE
table(CP) = CENT/VC;

'

mod <- mread(code=code_mod, model="ex2_mrg_poped")

data <- ev(amt=c(4,6,8,15), cmt=1) %>% as.data.frame %>% as.tbl %>% mutate(ID=1:4, dose = amt*10,amt=amt*10000)
out <- mod %>% data_set(data) %>% Req(CP) %>% carry.out(dose) %>% mrgsim(end=48, delta=0.1)

plot(out)

ff <- function(model_switch, xt, p, poped.db){
  times_xt <- drop(xt)
  dose_times <- 0
  time <- sort(unique(c(times_xt,dose_times)))
  is.dose <- time %in% dose_times

  data <-
    dplyr::data_frame(ID = 1,
                      time = time,
                      amt = ifelse(is.dose,p[["DOSE"]], 0),
                      cmt = ifelse(is.dose, 1, 0),
                      evid = cmt)

  out <- mrgsim(mod,end=-1, data=data,param=as.list(p),recsort=4)

  y <-  out$CP

  y <- y[match(times_xt,out$time)]

  return(list(y=matrix(y,ncol=1),poped.db=poped.db))
}

## Define other functions


## -- parameter definition function
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                VC=bpop[2]*exp(b[2]),
                Q=bpop[3],
                VP=bpop[4],
                TVLAG=bpop[5],
                DOSE=a[1],
                WT=a[2])
  return( parameters )
}
## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  y = y*(1+epsi[,1])+epsi[,2]
  return(list( y= y,poped.db =poped.db ))
}


poped.db <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175),
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1),
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=1,
                                  groupsize=6,
                                  xt=c(3,5,8,18),
                                  minxt=c(3,4,7,9),
                                  maxxt=c(4,7,9,24),
                                  bUseGrouped_xt=0,
                                  a = c(40000,70),
                                  maxa=c(40000,70),
                                  mina=c(40000,70))


system.time(p1 <- plot_model_prediction(poped.db))
p1


system.time(p2 <- plot_model_prediction(poped.db,IPRED=TRUE,IPRED.lines=FALSE, IPRED.lines.pctls=FALSE,
                                        DV=TRUE, DV.lines=FALSE, DV.points=FALSE, separate.groups=FALSE, PRED=TRUE,
                                        sample.times = TRUE, sample.times.IPRED = FALSE, sample.times.DV = FALSE))
p2


### Evaluate FIM
FIM <- evaluate.fim(poped.db)
FIM
det(FIM)
get_rse(FIM, poped.db)

## here we see that lots of parameters are poorly estimated
## not surprising, 1 group, 4 samples, and 4 bpop and 2 d parameters to estimate
poped_db_2 <- create.poped.database(ff_file="ff",
                                  fg_file="sfg",
                                  fError_file="feps",
                                  bpop=c(CL=20.9,VC=70.1,Q=1,VP=27.5,TVLAG=0.175),
                                  notfixed_bpop=c(1,1,1,1,0),
                                  d=c(CL=0.08,VC=0.1),
                                  sigma=c(0.06,6.48),
                                  notfixed_sigma=c(0,0),
                                  m=2,
                                  groupsize=3,
                                  xt=list(c(3,5,8,18),c(1,6,50,120)),
                                  minxt=0,
                                  maxxt=120,
                                  a = c(40000,70))
FIM <- evaluate.fim(poped_db_2)
FIM
det(FIM)
get_rse(FIM, poped_db_2)


# run only methods that can be parallelized to make things faster
poped_optim(poped_db_2, opt_xt = T,parallel = T, method=c("ARS"))  # approximate
#poped_optim(poped_db_2, opt_xt = T,parallel = T, method=c("ARS","GA","LS")) # more accurate

Hi, andrewhooker! In your code you do not use $MAIN double CL=TVCL*pow(WT/70,0.75)*exp(ECL); So how the WT work as covariate?

Change：OK, I know the answer. Here is an EX

ff <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    times_xt <- drop(xt)
    dose_times = seq(from=0,to=max(times_xt),by=TAU)
    time <- sort(unique(c(times_xt,dose_times)))
    is.dose <- time %in% dose_times
    data <- dplyr::tibble( ID = 1,
                           time = time,
                           amt = ifelse(is.dose,parameters[["DOSE"]], 0), 
                           cmt = ifelse(is.dose, 1, 0), 
                           evid = cmt,
                           CL = parameters[["CL"]], V2 = parameters[["V2"]], V3 = parameters[["V3"]],
                           V4 = parameters[["V4"]], KA = parameters[["KA"]], Q3 = parameters[["Q3"]]*(1+parameters[["Q3_ALT"]]*(parameters[["ALT"]]-15.2)), Q4 = parameters[["Q4"]],
                           F_DEPOT = parameters[["F_DEPOT"]], Q3_ALT = parameters[["Q3_ALT"]]
    )
    # times <- sort(c(times_xt,dose_times))
    out <- mrgsim(mod, end=-1, data=data, param=as.list(parameters), recsort=4)
    y <-  out$IPRED*(1000*F_DEPOT)
    y <- y[match(times_xt,out$time)]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

So, we can do like CL = p[["CL"]]*WT**CL_WT

from poped.