What are possible reasons for differences in numerical output from the same procedures performed in different statistical packages?

My colleague and I do the exactly same thing with the exactly same data and get different results. He uses Stata, I use R. We both use 32-bit Windows machines.

Are there any possible explanations apart from human failure for this?

The concrete problem is calculating c-indices for Cox models. His are up to 0.005 higher.

-
If you post an example of the actual functions or code that are resulting in differing results, we (or StackOverflow) might be able to provide more specific answers, although @ocram's answer is correct in general. – jthetzel Aug 9 '12 at 16:33

One possible reason might be default options. In particular, with Cox models, there are different methods for handling ties. Nearly all Cox regression packages use the Breslow method by default, but R uses the Efron method. Iteration limit and other control options might also differ from one package to another. You should also check whether STATA fits the proportional hazards model or the accelerated hazard model.

-
Although for my specific problem, that doesn't seem to be it, that's a very good point. – miura Aug 10 '12 at 7:02

At the CPU level, the result of adding a bunch of numbers up can depend on the order you do it. I had this happen to me when doing some OpenMP computations. The result wasn't constant, even though there was no randomness in the program. What was happening was that the order in which the OpenMP threads finished was variable, and so the sum computed was being done in a different order, and giving a different result. I eventually tracked this down to a minimal C example that exhibited the behaviour - add the vector up from the first to the last and you got a different answer to if you added it up from the last to the first.

Although the differences were tiny, probably one part in 10^10, it had a large effect on my eventual output. This was a maximisation problem, and the tiny difference in slope was sending the maximiser to a different point in the space. After 100 iterations, the maximiser could well be in a very different place just because of the different ordering of an addition in one case. And once the system has moved, there's no way back (especially when the likelihood surface is flat...). Its a sensitive dependence on initial conditions (chaotic) situation.

So, even if your programs are implementing the same algorithm correctly, they may still get two different answers. Both are right in some sense, both are wrong in some sense.

Whether they are significantly different is what's important. Don't sweat the 0.00001 of a probability difference...

Of course with R you can look at the source code and find out how it does it...

-
I suspected something like this, but the differences I observe are much larger than 10^-10. They are in the magnitude of 10^-3, which is considerable for a statistic that is bound between 0 and 1. – miura Aug 10 '12 at 7:03
The issue I had was that although the difference in computation of the likelihood was only 1e-10, it had a substantial difference in the final maximised parameter values - 26.1, or 45.7, or -123.5. Small differences in part of an underlying algorithm can have chaotic effects. – Spacedman Aug 10 '12 at 8:26