Sigma is 0 when calculating hess_fnc

[elena-krismer]

Hi, Our prottipackage is using proDAfor the differential abundance calculation. One of our users recently encountered following error, when using protti:

proDA_fit <- proDA(proDA_input_clashing, design = proDA_design_clashing)

After further investigation I found that:

• When fitting the probabilistic dropout model, the algorithm ensures that sigma never falls below 1e-100, see In line 280 pd_lm.R - however R determines 1e-100284 as the bigger number (for what reason ever??) and assigns it to sigma2
• When passing sigma2 to the hess_fnc 1e-100284 gets lost on the way and sigma2 will be read as 0 .
• This ultimately results in a division by zero, leading to ‘ fit_sigma2_var’ in line 322 of pd_lm.R becoming ‘NA’
• Consequently, the check if( fit_sigma2_var < 0) fails

I guess the error comes from the fact that R is very bad when it comes to handling values very close to zero and can only handle doubles/64 bit floating numbers. When comparing 1e-100284 (in this case) with 1-e-100, R messes up.

[elena-krismer]

here would be the link to the data to reproduce the error: https://github.com/elena-krismer/proDA_troubleshoot_data

[const-ae]

Hi Elena,

this is an excellent bug report; thanks a lot for taking the time to investigate what went wrong and how it can be fixed and making a minimal PR.

Checking for sigma2 <= 0 in hess_fnc seems reasonable to me. Is there a particular reason though that it is preferable over changing line 314 to fit_sigma2 <- max(1e-100, nl_res$par[p+1])? This would mirror the behavior in lines 280, 290, and 300.

Alternatively, we could keep your check in hess_fnc and remove the checks in lines 280, 290, and 300. What do you think is best?

Best, Constantin

[const-ae]

Thanks again for your contribution. I merged the PR.