How can I calculate the Hessian matrix for fit results with sbiofit (SimBiology)?

Hi all,

With sbiofit (SimBiology) the fit results give me standard errors obtained from the Jacobian calculation. I would like to compute the standard errors with the Hessiana matrix, but I don't know how to get it analytically from my object function and the fit parameters.

Thank you in advance,

Valentina

Answer 1

[Arthur Goldsipe]

I can't think of an easy way to calculate the Hessian matrix. The objective function is not analytical (since it involves approximate numerical solutions of differential equations). It's possible to estimate the Hessian numerically for situations like this, but I worry that it would be difficult to get an accurate estimate. (We already have challenges calculating the Jacobian accurately.)

Have you looked at the different options available for calculating confidence intervals using sbioparameterci or sbiopredictionci? If none of those are satisfactory, can you tell me more about what you're looking for? Can you point me to any references that describe the calculations you're trying to do? I'd like record this in our enhancment request database, and the more details the better.

-Arthur

[Valentina Vasic]

Yes I realized that the calculation with the Jacobian is not accurate, but my problem is different. I also tested these two functions and the confidence intervals give me more or less the same standard errors calculated difectly by sbiofit (but not equal). About these two functions. I wanted to know if they take into account the standard errors of the output of sbiofit or do they calculate the confidence intervals from zero without considering the standard errors of sbiofit?

Getting back to my problem with Hessian matrix. I am using fits with sbiofit and trying to fit my data also with Bayesian information. I have modified the object function and the values I get from the fit of my parameters are compatible with values found with another software where the use of Bayes values is already implemented in the search for minima. The estimated parameters therefore turn out to be the same between the two software, but the standard errors are not. With sbiofit I get larger errors than with the starting software.

Following the calculation of standard errors, these are done through the Jacobian (line 1826 in SimBiology.fit.internal.FitObject version of Matlab 2021a). Unfortunately, Bayesian parameters are not taken into account in the Jacobian calculus. This causes the values to be small, which then inverted give rise to large standard errors. In particular, for one of my cases, I have 15 measured points and 17 parameters to estimate. Without the Bayes, I cannot estimate, but knowing the Bayes of these 17 parameters, I should get good results. This given that I consider the Bayes as additional measured points. So if the Jacobian does not consider the Bayes, I get standard errors of the order of 10^4/10^5. The idea is to use the Hessiana matrix and calculate it on the object function with respect to the parameters to be estimated. Avoiding the Jacobian calculation and finding the errors by inverting the Hessian matrix. Since it is calculated on the objective function, it contains the Bayesian information and I should solve my problem.

Now the problem is how to insert the Hessiana calculus in row 1826 in SimBiology.fit.internal.FitObject (Matlab 2021a)? Or however to compute the Hessiana to correct the computation of covBeta, sse and mse (line 1769 in SimBiology.fit.internal.FitObject (Matlab 2021a)) and consequently the standard errors. Basically, I wanted to call the object function from line 293 in SimBiology.fit.internal.NegLogLikelihood as the model is called in line 1821 of SimBiology.fit.internal.FitObject.

I hope I have provided you with more information and that my problem is clearer. Thank you so much for your help.

[Arthur Goldsipe]

I'm not sure I completely understand all of the details. I don't understand how you can accurately estimate the Hessian matrix if your Jacobian is inaccurate.

Are you saying that you modified the SimBiology code to customize your objective function? If so, that could explain the incorrect Jacobians you are seeing. By default sbiofit will attempt to calculate the Jacobian rather than relying on the optimization function to estimate the Jacobian with finite differencing. The Jacobian calculations that sbiofit does would need to be updated if you make any changes to the objective function calculations.

If disabling sensitivity analysis gives you accurate Jacobians, then one idea I can suggest is to use fmincon for parameter estimation, since it reports an estimate of the Hessian matrix in its results. You can find more info here.

I'd also like to understand more about your Bayesian approach. Can you point me to any references on this approach? It sounds like this would be useful functionality to add to SimBiology.

[Valentina Vasic]

Of course the Jacobiana is wrong, since the function it takes as imput is the model itself, but not the object function. If I could pass the object function to it, instead of the model, then the calculations would be correct.

The point is that I do not know how to pass the object function to it.

From theory I found a method to estimate the covariance matrix from the Hessiana (calculated with the object function). This should work and be more accurate than the calculations made with the Jacobian. The Hessiana should obviously be calculated numerically, similar to the Jacobiana.

In particular, what would you like to know about the Bayes method? I have followed some information from this book.

[Arthur Goldsipe]

I'm sorry, but I still don't fully understand what you're trying to do. I think we're at the point where I need to see sample code before I can recommend anything more specific. If you're interested in that, please feel free to share code here or to reach out to me via my profile if you want to share code privately.

I'm familiar with Numerical Recipes in C, but I don't recall anything in it about Bayesian approaches. What I was hoping to get was some sort of reference that describes whatever algorithm or approach you're trying to implement, in sufficient detail that I could think about how to implement it in SimBiology.

[Valentina Vasic]

I have forwarded you an email privately. And I would be happy if we could talk by email.