kylebaron / klava

Parameter optimization with mrgsolve
5 stars 0 forks source link

ELS method #7

Open Mutaz94 opened 4 years ago

Mutaz94 commented 4 years ago

When trying to fit individual data using els method, I face this problem. Any particular suggestions?

library(dplyr)
library(mrgsolve)
library(nloptr)
library(ggplot2)
library(rlang)
library(klava)
data <- read_csv("data.csv")

mod <- mread("model1.cpp")

data <- data %>% filter(CMT %in% c(1, 2)) %>% mutate(AMT = replace(AMT, AMT==".", 0)) 
data <- data %>% 
  mutate(EVID = if_else(CMT==2, 0, 1))

ggplot(data, aes(TIME, DV))+ geom_point()+theme_bw()

theta <- all_ident(tvcl = 18, tvv = 30, tvkatr = 6.5,
                   tvic50 = 0.22, 
                   tvcort_in = 1.66,
                   a0 = 0.44,
                   a1 = 0.22,
                   b1 = 0.16,
                   a2 = 0.04,
                   b2 = -0.05,
                   pi = 3.14159265359,
                   hill = 3)
data <- data %>% mutate(AMT = as.numeric(AMT))
data <- data %>% mutate(CMT = if_else(CMT==2, 0, 1))
fit <- fit_nl(theta, data, mod = mod, pred_name= "Conc", cov_step=F,
              pred_initial=TRUE, ofv = els)

Error message:

Fitting with els ...Error in res^2/sigma : non-numeric argument to binary operator

the data file contains: ID, TIME, DV, CMT, MDV, WT

for the model file


[init]
gut = 0 ,
  cent = 0.44,
  endo =0.44,
  trans1 = 0,
  trans2 = 0,
  trans3 =0

[param] 

tvcl =22.6 
tvv = 41.2 
tvkatr = 10.5
a0 = 0.44
a1 = 0.22
b1 = 0.16
a2 = 0.04
b2 = -0.05
pi = 3.14159265359
hill = 3
tvic50 = 1.02 
tvcort_in = 1.66

wt = 20 

[omega]
0.0818 0.0431 0.15 0.111 0.941 
[sigma]  
0.0326 0.0098 
[main] 

double cl = tvcl*pow((wt/70), 0.75);//*exp(ETA(1)); 
double v = tvv*(wt/70);//*exp(ETA(2)); 
double kel = cl/v; 
double katr = tvkatr;//*exp(ETA(3));
double ic50 = tvic50;//*exp(ETA(4));
double cort_in = tvcort_in*pow((wt/70), 0.75 );//* exp(ETA(5));

double s2 = v/100;

[ode] 

#define t SOLVERTIME

double rcort = a0*kel+((a1*kel+b1*(2*pi*1/24))*cos(2*pi*t/24) + (b1*kel+a1*(2*pi/24))*sin(2*pi*t/24) + (a2*kel+b2*(2*pi*1/12))*cos(2*pi*t/12) + (b2*kel+a2*(2*pi/12))*sin(2*pi*t/12));
if(rcort<=0) rcort = 0.001; 
double cp = cent/s2 ;
double effect = pow(cp,hill)/(pow(ic50,hill)+pow(cp,hill)); 
//double effohp = pow(cp,hill)/(pow(i50,hill)+pow(cp,hill));
//double effd4a = pow(cp,hill)/(pow(d50,hill)+pow(cp,hill));

//------------ Differential Equations-----------------//

dxdt_gut = -katr*gut;
dxdt_trans1 = katr*(gut-trans1);
dxdt_trans2 = katr*(trans1-trans2);
dxdt_trans3 = katr*(trans2-trans3);
dxdt_cent = katr*trans3 - kel*cent; 
dxdt_endo = cort_in*rcort*(1-effect)-kel*endo; 
//dxdt_ohp  = rohp*rcort*(1-effohp)-kout_ohp*ohp;
//dxdt_d4a  = rd4a*rcort*(1-effd4a)-kout_d4a*d4a;

[table] 
double Conc = (cp + endo)* (1+EPS(1)) + EPS(2); 
//double OHP = ohp;//*(1+EPS(3)); 
//if(OHP <= 0) OHP = 0.001;
//double D4a = d4a;//*(1+EPS(4)); 
//if(D4a <= 0) D4a = 0.001; 
if(Conc <= 0) Conc=0.001;
[capture] 
Conc 
kylebaron commented 4 years ago

I can't run the example, so sort of hard to tell. Can you try adding a "sigma" parameter to theta? Also, are you estimating the value of pi there?

Mutaz94 commented 4 years ago

Thanks! The algorithm worked fine when I insert the sigma in the thetas. if I want to insert 2 sigmas, should we call it separately or like two sigma section. how mrgsolve will look at it in combined error model?