why PSI values are diffrent when nlmefit using ode45 solver compared to when it uses analytical solution?

1 view (last 30 days)
Sepi Sharif on 8 Sep 2021
Commented: Sepi Sharif on 20 Sep 2021
nlmefit estimated parameters differently when it uses analytical solution compared to when it uses ode45 solver for the same dataset. PSI values are completely underestimated.
[beta,PSI,stats,b] = nlmefit(timeNew,DVnew,IDnew,dosePerID,model,logparameter0,...
'REParamsSelect',[1 2 3 ],'ErrorModel','Constant','ApproximationType','FOCE','CovPattern',P,'RefineBeta0','on',...
'ParamTransform',repelem(1,3),'OptimFun','fminunc');
function concentration=Getsolutionode(parameter,t,doseNew)
%% initial values
V=parameter(:,1);
ka=parameter(:,2);
CL=parameter(:,3);
predictor(:,1)=t;
minTime=min(t);
maxTime=max(t);
y0=[doseNew 0];
%% ode45
options=odeset('RelTol',10-9,'AbsTol',10-9);
solode=ode45(@odefunction,[minTime maxTime], [y0],options,V,ka,CL);
rspl=deval(solode,t);
concentration=rspl(2,:)'./V;
end
function dxdy = odefunction(~,y,V,ka,CL)
kel=CL/V;
dxdy= [-ka*y(1)
ka*y(1)-kel*y(2)];
end
Sepi Sharif on 20 Sep 2021
Hi all,
a question, about nlmefit, approximation type: so, dose LME as the default mean that nlmefit also, can solve linear mixed effects as well? say a systme of linear ODEs could be solved by nlmefit?

Star Strider on 8 Sep 2021
From the description, it appears that the code is attempting to fit data to the solution of a differential equation, or a system of differential equations.
To do this with ode45 (or more appropriately ode15s for kinetic equations), see Coefficient estimation for a system of coupled ODEs for a useful illustration of a similar problem.
.
Star Strider on 8 Sep 2021
As always, my pleasure!
I see that Arthur Goldsipe has responded, so I will leave you in his capable hands.
.

John D'Errico on 8 Sep 2021
Edited: John D'Errico on 8 Sep 2021
There can sometimes be problems when you try to apply a nonlinear estimation scheme on top of an ODE solver. The nonlinear fitting tool presumes that the objective function is a differentiable function of the parameters. But is it? Remember that an ODE solver is a tool designed to work with a tolerance. Change the parameters slightly, and you will get a subtly different result, but only to within a tolerance. Effectively, your objective function is not a differentiable function of the parameters, because there is an arbitrary and unpredictable fuzz on top of your objective.
How is this different when you use an analytical solution to the problem? Now the objective will generally be a smooth, differentiable function of the parameters.
Effectively, this CAN be a problem anytime you layer one numerical solver on top of another. One thing to try is to crank down on the tolerance of the internal solver. But will that always be sufficient? Possibly not, since if the problem is a stiff one, then the solution can vary significantly as you cary the parameters.
Another trick may be to use a stiff solver. ODE45 is a tool that is not designed to solve stiff problems. So you may gain by simply moving to a stiff solver. That is, solvers where the name ends in s, like ODE15s, or ODE23s.
Can I know that this is the reason? Not without a concerted effort to dig into exactly what you did. Is it possible that your ODE may be a little stiff? Possibly. That would certainly exaggerate any issues in the solve. And of course, there may be a bug in the code you wrote. :) Surely not! That would never happen. ;)
Sepi Sharif on 8 Sep 2021
analytical solution with nlmefit has these estimated parameters
V=8.536; ka=0.585; CL=0.140, PSI=[0.062 0 0
0 0.348 0
0 0 0.083 ]
But;
ode45 solution has
V=8.827; ka=0.772; CL=0.171 , PSI=[0.048 0 0
0 0.0000 0
0 0 0.0001]
there is a discrepancy here!
So, the analytical solution is
concentration=doseNew.*ka./(V.*(ka-kel)).*(exp(-kel.*(t))-exp(-ka.*(t)));
where, kel=CL/V;
this is a nonlinear_nonstiff problem!

Arthur Goldsipe on 8 Sep 2021
The reason the results from nlmefit differ significantly is because your ODE solution is not equivalent to your analyitical solution. I think you need to update file Getsolutionode.m to set minTime to 0 rather than min(t). Otherwise, your initial dose condition gets applied at min(t) instead of at time 0, and the ODE solution essentially has a time-delay when compared to your analytical solution.
In addition, as I mention in my comments on the question, I recommend changing 10-9 to 1e-9 when setting the tolerances. Otherwise, you will be using extremely loose tolerances when solving the ODEs.
Finally, I see from your example code that you are looking at a standard problem from pharmacokinetics (the "warfarin" dataset). This is exactly the sort of problem that SimBiology is designed to solve. You can build your pharmacokinetic model and run nlmefit and other analyses using SimBiology's graphical user interface or using command line functions and scripts. Here are some advantages I see to using SimBiology to solve this sort of problem:
• You can build your models graphically.
• You can select many standard PK models (including this one) from a built-in library.
• You can represent your model in a more natural way, using compartments, species, and reactions.
• You don't have to worry about the "bookkeeping" required to write the underlying differential equations.
• You can use built-in programs/tasks in the graphical user interface to perform mixed effects modeling.
• Your simulations almost always run faster. SimBiology can solve linear ODE models like this one analytically, using highly optimzed C++ code. And models that can't be solved anlytically can be solved much faster by taking advatange of its acceleration feature.
If you have a SimBiology license, I can show you how I would solve this problem in SimBiology. Just let me know whether you prefer to use the graphical user interface, the command line, or both.

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!