-
Notifications
You must be signed in to change notification settings - Fork 20
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
troubleshooting tips - Zero FIM
determinant
#2
Comments
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,
|
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 :/ |
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
fg and feps
poped.db
|
still trying to diagnose and going through the functions one at a time. I see that for some reason the UPDATE: never mind - realized that |
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 ( Can you post the model code too? |
Thanks, Kyle - here is the model code.
|
First: Second: If it was just the first issue, we could take this off the If you have:
in the data base, does |
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
|
Thanks @kylebmetrum. A couple of changes. You had |
@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 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 |
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 |
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. |
Hi, andrewhooker! In your code you do not use Change:OK, I know the answer. Here is an EX
So, we can do like |
hi Andy,
I am not sure how to diagnose the condition where the
evaluate.fim
returns a0
determinant. I know that myff_file
returns a validy
vector (checked by running the setup outside a function setting). Thepoped.db
seems to be have been set up correctly as far as the design goes, this I can confirm by plotting the predictions using theplot_model_prediction
which returns a population prediction at the specified design times. But when I run anevaluate.fim
on thispoped.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 thatA
andB
below are parameters that KA is derived from.The text was updated successfully, but these errors were encountered: