I am unsure if I am misreading the documentation, if this is a bug in the documentation or a bug in the code. In any case the documentation for nfbasis suggests that if I call:
nfbasis([x^4 - x^3 + 10558*x^2 - 2169592*x + 195873967, [13,73,89]])
I should obtain a global integral basis because:
"The result is actually a global integral basis if all prime divisors of the field discriminant are included!"
However, even though the discriminant of the field in this case is (13*73*89)^3, the basis returned has index (103*499), so is not globally maximal.
For reference, I am using gp version 2.11.1.