`fit_alpha_mle` questions

[mikelove]

in the loss and dloss functions, why not just use log_alpha here:

https://github.com/owkin/PyDESeq2/blob/main/pydeseq2/utils.py#L612

and here:

https://github.com/owkin/PyDESeq2/blob/main/pydeseq2/utils.py#L633

You can likewise just define log_alpha_hat once to avoid recomputing it on every function and gradient evaluation.

And since this is optimizing in the log space, I wouldn't put an alpha in the denominator (we just differentiate wrt log_alpha). Instead of:

reg_grad += (log_alpha - log_alpha_hat) / (alpha * prior_disp_var)

I would have:

reg_grad += (log_alpha - log_alpha_hat) / prior_disp_var

I'm not sure if this is why, but I see a lack of shrinkage sometimes for high counts (where we do shrink these in DESeq2), and I'm trying to figure out the source. Have you seen this?

[BorisMuzellec]

Hi, indeed there are some inconsistencies in the optimization. It was designed to optimize w.r.t. log(alpha) but the derivative was computed w.r.t. alpha. I opened a PR (#107) to fix this.

You're absolutely right about unnecessary recomputation of log(alpha) and log(alpha_hat). Computation of W could also be avoided when CR regularization is not applied.

Do you have a link to some code (or a piece of code) to reproduce the figure above, so I can see whether it fixes the issue?

Thanks again for your help in improving the code!

[mikelove]

Impressive work implementing everything.

I do have code to reproduce something like the figure here:

https://github.com/mikelove/test_pydeseq2/blob/main/test_pydeseq2.R#L7-L11

The tricky thing is, I'm not sure if this is the cause, as you have the correct function definition, but I think incorrect gradient (given what we are taking derivative of is log alpha not alpha). Depending on the optimization method, it can get around an incorrect gradient definition. But maybe this clears things up. I will also test on my end.

[mikelove]

This appears to be fixed via #107

[BorisMuzellec]

Thanks a lot for checking that #107 fixes the problem!

I'm closing this issue