Yes, but as long as the imaginary part is there, it is impossible to rule out a doubt (unless you performed other converging methods independently, which you did :-)
As far as the Mellin/InverseMellin transform method is concerned, the convincing workaround would rather be:
InverseMellinTransform[-Tan[ Pi \[Rho]/2 ] MellinTransform[
Exp[- t^4], t, \[Rho]], \[Rho], x] // FullSimplify
(See Frederick W. King, Hilbert Transforms, Cambridge U.P, p. 274, which works for even functions).
This yields the elegant :
(1/2 - I/2) E^-x^4 ((-1 + I) + GammaRegularized[1/4, -x^4] -
I GammaRegularized[3/4, -x^4])
Now for the reduction to reals, I have not completed the algebra yet, but the reasoning goes as follows: incomplete Gammas with a strictly negative argument are rotated in the complex plane by Pi theta / 2 with respect to the value for positive arguments (see NIST handbook, 8.5.5, Confluent Hypergeometric Representations), where theta is the form index.
Here the increase in the form index is 1/2 for an argument that can only be strictly negative or null, hence a rotation by Pi /4 from the first Gamma to the second. The null case is obvious, and otherwise after some basic algebra it follows that the big bracket is proportional to Exp[ i Pi / 4], which multiplied by the (1 - I) initial factor yields a real.
