Standard errors of parameters fitted with sbiofit using Bayes are too high

Hi,

I am currently working with SimBiology and have implemented a compartmental model to follow the Lutatera therapy. I am using the 2021 version of Matlab. I have modified sbiofit to implement the object function in order to perform a fit with Bayes values. Comparing it with another software that already uses Bayes, I could see that the fitted parameters have values that differ little. The problem is that with matlab, the standard errors that sbiofit calculates are too large. At first I thought it was just because of the degrees of freedom, but unfortunately, even when correcting them, I still get too high errors on the fitted parameters. Upon redoing the code, I noticed that the Jacobian does not take into account my Bayes values, nor the standard Bayes error. numjac is the function that calculates the Jacobian.

My question is whether it is possible, and how, to incorporate Bayes standard errors into the standard errors of my parameters. In theory, from the fit, the errors of the fitted parameters should be no larger than the Bayes errors. Has anyone ever solved this problem in other functions used for the fit? Or are there functions that already use Bayes that can help me solve my problem?

Regards

[Florian Augustin]

Hi Valentina,

SimBiology's function sbiofit currently (R2022b) does not support Bayesian approaches for parameter estimation. If you would like to implement Bayesian parameter estimation, you could take a look at mhsampler or hmcSampler in the Statistics and Machine Learning Toolbox. I am planning on making a FileExchange submission for Bayesian parameter estimation using mhsampler with SimBiology models soon. I can post a link here if you are interested.

To your actual question. Can you post more details about the algorithm you implemented? Any reference, code that we can run to reproduce the issue would help the community to answer your question.

Best,

Florian

[Valentina Vasic]

Dear Florian,

thank you very much for your reply, I will have a look at these two functions.

I don't want to post the whole code, because it is too big. Basically in the NegLogLikelihood.m I modified the object function by adding Bayes values as additional measured points, weighted. At the time I had already found questions posted about this and a colleague answered them in great detail. So the fitting parameters are correct for me.

What is wrong are their standard deviations (SD) of them. Their calculation comes at the end of the fitting in the function FitObject.m. At first I got SD = Inf or SD=NaN. That was due to the degrees of freedom that were calculated dfe = max(0, (numValidResiduals-rankJ)) in the function calcCovBeta ([sse, mse, covBeta] = calcCovBeta(jacobianStacked, residualStacked, beta)) of FitObject.m. The function is located towards the end of the code.

Now, when I estimate the parameters-SD, there are relatively high compared to the SD that other software gives me and works in similar way. I thought that in the calculation of the Jacobian (getjacobian function in FitObject.m) the Bayes values must also be entered, but I don't know how. If my Bayes are as simple measured points added to the problem, in theory the Jacobian should have a dimension of (N+B)x(N+B), where N are the number of measured points and B the number of Bayes. But the dimension is NxN. So to solve my problem I think I need to make the Jacobiana larger or as large as I can make it take Bayes values into account.

Best,

Valentina