I was looking to run Tukey HSD with statsmodels, and when I compared the results to R, I was quite startled to see very different numbers. Digging a little further, I noticed that the difference is from the statsmodels.stats.libqsturng approximation to the distribution.
I found some pascal code and ported it to python
Any thoughts if/how to proceed?
Sorry but newbie contributor here. I put the code in my fork, which I didn't realize would suggest creating the PR in the main repo.
Other differences with statsmodels.stats.libqsturng:
* Pascal code authors validated the correctness of their code in the paper, and is not an approximation
* libqsturng only supports p-values from 0.1 to 0.999
* libqsturng is same as R code _only_ up to 1 or 2 digits
I believe the difference is because libqsturng uses a table of pre-calculated values,
using the same algorithm as R, but the table has different values than R.
Since I don't know the code or the paper for libqsturng, I hesitate to update that table
which might have unexpected results.