Initial "seed" for maximum likelihood estimation

[Tymon Sloczynski]

I have a fairly complex do-file that estimates multiple binary choice models using maximum likelihood (especially probit, using the probit command). In a small number of those models, I see minor differences in estimation output depending on whether I run this code immediately after launching Stata or after running some other code. Only as a result of this did I learn that maximum likelihood estimation in Stata is not always exactly reproducible when the same code is run repeatedly.

I want my code to deliver the same output as it does immediately after launching Stata, regardless of when it is run. Is there a special "MLE seed" that would accomplish this?

For example, when typing

set seed 123456789

I can replicate the pseudo-random numbers Stata generates immediately after launching, even though my current session has lasted much longer.

Similarly, when I type

set sortrngstate 1001XZA112210f4b16c1cb10507a1f38cb440c40003c9a8356 6fa1201b69ab0aff2f2401e18

I can replicate the sort order obtained immediately after launching Stata.

Is there a similar "seed" that ensures maximum likelihood estimation runs as if Stata had just been launched?

[Joseph Coveney]

I have a fairly complex do-file that estimates multiple binary choice models using maximum likelihood (especially probit, using the probit command). . . . I see minor differences in estimation output depending on whether I run this code immediately after launching Stata or after running some other code. . . . maximum likelihood estimation in Stata is not always exactly reproducible when the same code is run repeatedly.

To my knowledge the only conditions under which Stata's maximum likelihood estimators would have any random behavior is when the default initial values give rise to a missing-valued log-likelihood and a random search is undertaken in order to find feasible starting values.

But even then, the likelihood function for a conventional (single-level / nonhierarchical) binary regression model is strictly concave and so it shouldn't matter where the random starting values begin; they should all migrate to the same global maximum. (You can always double-check that by profiling the likelihood and seeing whether it's multimodal or otherwise kinky, but if that's the case—complete or quasicomplete separation, collinearity—then the corrective action requires something other than setting the seed.)

So, my take on your situation is that you have some indeterminacy in your "fairly complex do-file", most likely in a data-munging step involving a sort or bysort command. Here, too, you'd be better off by making the data-management process inherently determinant.

Is there a similar "seed" that ensures maximum likelihood estimation runs as if Stata had just been launched?

Why would the state of the random-number generator have to be as if Stata had just been launched? Wouldn't if suffice if it is just the same each time that the do-file is run? For that, you can choose any arbitrary seed value and institute it in a set seed command placed in the beginning section of your do-file.

I suppose that if you really had a reason for insisting on ensuring that the random-number generator's state is exactly that for when Stata is just launched, then you could open a fresh instance of Stata and

file open file_handle using initial_state.txt, write
file write file_handle "`c(rngstate)'"
file close file_handle

and then at the top of your do-file insert the following

file open file_handle using initial_state.txt, read
file read file_handle state
file close file_handle
set rngstate `state'

But, again, I recommend looking for the source of your reproducibility problem elsewhere, first in the data-management code (convert all sort . . . commands to isid . . ., sort etc.) and then in the profile likelihood of the problematic probit regression models.

[Tymon Sloczynski]

Thanks a lot! At this point, it appears that setting sortrngstate to its startup value resolves the issue. I need to analyze this further, but at this point, I am optimistic that this is the solution.

[Joseph Coveney]

. . . it appears that setting sortrngstate to its startup value resolves the issue. . . . at this point, I am optimistic that this is the solution.

You might have misunderstood me. Setting the sort seed is the one action that I specifically did not mention—it does not make data-management process inherently determinate.

To quote Clyde Schechter:

It is true that irreproducibility is often the result of indeterminate sorting. But -probit- should not be subject to that problem: the results really should be the same regardless of the sort order of the data. It is more likely that something you are calculating before you get to the -probit- command depends on the sort order. The solution, then, is to fix that calculation so that it is independent of the sort order, not to stabilize the sort. Stabilizing the sort can be done easily, but it just sweeps a huge problem under the rug.

There is even more from StataCorp itself:

Stata/SE and Stata/MP use the jumbler differently and so produce different orderings of ties, even when starting from the same seed/state. What’s more, Stata/MP with two processors and Stata/MP with four processors also produce different orderings of ties. The older qsort (prior to Stata 17) and the newer fsort (Stata 17 and beyond) also use the jumbler differently and produce different orderings of ties, even when starting from the same seed/state. . . . So any reproducibility produced by set sortrngstate is specific to the edition of Stata that you are running and which sort method is being used.

The approach you're taking risks your not being able to reproduce your own analysis's results even within Stata in the future. Now can you imagine what someone in your audience who is using R (or SAS or whatever) to reproduce your analysis will conclude?

I recommend making the data-management process inherently determinate, that is, not dependent upon particular sort order that obtains with a given sortrngstate in a given flavor and release of a particular software package.

[daniel klein]

It should never be necessary to set the sort seed. If results change because Stata randomly breaks ties, that is strong evidence that some relevant information has been omitted from the sort key. If all relevant sorting information were included, the ordering would be fully determined and the results would be stable. Setting the sort seed does not solve the underlying problem; it merely masks it.

[Sebastian Kripfganz]

That's interesting, in a couple of ways.

As Joseph mentioned, randomness with maximum-likelihood estimation would usually come from different initial values. While the estimator should converge to the same solution if the likelihood function is well-behaved, numerical results might still differ depending on the chosen precision for the convergence criteria.

If the results are substantially different, then this would suggest an ill-behaved likelihood function with multiple local maxima or flat regions. In this case, fixing a seed would make the results reproducible, but this would actually only mask the fundamental problem. But I am sure Tymon Sloczynski knows all of that.

Different sort orders should normally also not matter for the final result. If they do, this again suggests that there might be a more fundamental issue with the code, and fixing the seed might just hide this problem.

What really surprised me, though, was that Stata indeed appears to have an initial startup seed of 123456789. (It's documented in help set seed.) In effect, this makes random numbers produced by Stata "non-random". While they are only pseudo-random anyway, with this fixed startup seed, they are no longer "as good as random". I don't think Bill Gould would have approved that. I always was under the impression that the initial seed would be a function of startup conditions; e.g., the time at which the Stata instance was started. Apparently, that's not the case. Does anyone know a reason why Stata uses this fixed startup seed?

[Joseph Coveney]

. . . numerical results might still differ depending on the chosen precision for the convergence criteria.

Yeah, I did consider that before responding, but the default convergence tolerances for ml are reasonably tight, and wouldn't make a night-or-day difference in the regression coefficients.

Does anyone know a reason why Stata uses this fixed startup seed?

In the interests of reproducibility for those who neglect to set the seed, maybe: there's really no advantage to an Excel style RANDOMIZE TIMER approach under those circumstances.

[Tymon Sloczynski]

Thank you, all. I agree 100% that the seemingly random behavior of MLE is a major cause for concern; no need to convince me! However, in my specific context, this occurs in only a handful of regressions in a sequence of about one thousand, and I absolutely do not care about those particular regressions per se; I only care about the average behavior of a statistic I compute repeatedly across those many instances. What this means is that I am "sweeping under the rug" the fact that, in one run, I may conclude that something happens 74.9% of the time, and, in another run, it happens 75.3% of the time, while the only thing that matters for my conclusions is that it happens "pretty often." However, I also need to be able to reproduce the exact value I already reported, even if it doesn't matter in any substantive sense whether it's the one or the other.

[Joao Santos Silva]

Dear Tymon Sloczynski et al.,

This is may be unrelated to your problem, but one thing I noted is that certain estimation commands change the sort seed; it took me ages to find this out and fix the problem. The following example illustrates this

cls
sysuse auto, clear
query sortseed
di r(sortseed)
qui probit foreign price mpg i.rep78
query sortseed
di r(sortseed)

As far as I can tell, this only happens if there is separation/perfect predictors.

Best wishes,

Joao

[Tymon Sloczynski]

I actually think that this is likely very related to my problem. Thanks for a good example, Joao!