Scott Cunningham writes:

You’ll appreciate this I think. I ran Claude code on 96 specs for a popular difference-in-difference estimator with the identical specifications, ranging covariates only, for three languages (R, Python and Stata) and 2 packages for each. Different estimates — sometimes five times off. Appears to be coming from how each one handles failures on logits estimation with perfect separation, with some just refusing and others refusing and passing on garbage and others dropping down to run a logit on a constant despite being told to use Xs. And it has to do with how each package handles “big numbers” and this machine epsilon business that I barely know anything about. The fact that economists have pointed on selection on observables means they also haven’t invested in basic diagnostics around estimating propensity scores. But even turning Xs into z-scores still didn’t fix it. It’s the most extreme version of the garden of forking paths I think because people just assume languages and packages will all agree.

Ahhhh, perfect separation in logistic regression . . . one of my favorite topics! For an updated take, see section 14.6 of Regression and Other Stories. It makes sense that different packages will handle this differently.

Here’s Scott’s summary from his above-linked post:

The data come from a staggered rollout of mental health centers across Brazilian municipalities, 2002-2016. I have 29 covariates covering demographics, economics, health infrastructure, and population trends. Sixteen specifications sample the space from zero covariates up through all 29. All variation is coming from covariates, both within a package (i.e., which covariates) but more importantly across packages. And that’s the headline — even for identical specifications, you can get variation in estimates and sometimes as much as anywhere from an ATT of zero to an ATT of 2.38. Even if you drop the zero, some specifications yielded estimates ranging from 0.45 additional homicides to 2.38 homicides (measured as a municipality population homicide rate).

A five fold difference in estimated effects of a mental health hospital changing in the homicide rate is not trivial.

I guess one message is that there’s only so much you can learn from this sort of observational dataset. Ultimately it’s one data point in a larger explicit or implicit meta-analysis. More specifically, I like that Scott starts with the basic comparison and moves from there. I always think it’s a mistake to try to jump to the end.

Also I’d like to see a scatterplot with one dot for each municipality-year, plotting the outcome vs. some pre-treatment level, with different colors corresponding to treatment or control. I know that’s not the point of Scott’s post, but statistical packages can do graphs too! And more graphs, I’m not sure exactly what, but some way of getting a sense of the data. Ultimately this is being framed as a treatment vs. control comparison, so I want to see what the treatment and control data look like, before and after.

The multiverse!

Relatedly, Jessica points to two recent papers proposing the use of automated data analysis programs to perform multiverse analysis:

Many AI Analysts, One Dataset: Navigating the Agentic Data Science Multiverse, by Martin Bertran, Riccardo Fogliato, and Zhiwei Steven Wu

Editorial Screening when Science is Cheap, by Nic Fishman and Gabriel Sekeres

Jessica has discussed some related ideas at the end of this post and also here. And see here for lots more.

I agree with Scott and Jessica that these ideas all fit together somehow.