Category Archives: Multilevel Modeling

“Death and Taxes: Taking Causal Inference into the Real World”

Posted on by
1

Nandita Mitra is giving a webinar this Tuesday at noon with the above title and the following abstract:

Statisticians develop methods with the hope that they will ultimately inform decisions and improve the world around us. But taking causal inference into the real world requires more than choosing the right estimator: it requires identifying the question that matters, engaging deeply with its context, and communicating statistical reasoning to very different audiences.

In this talk, I will draw on two experiences in which causal inference took me well beyond traditional statistical applications. I will first discuss my experience serving as an expert witness in a death penalty case involving racial disparities in jury selection, where ideas including confounding, mediation, and individual- versus population-level effects had to be translated into compelling and understandable arguments in a court of law. I will then describe my involvement in evaluating Philadelphia’s excise tax on sweetened beverages, where cross-border purchasing led us to confront interference and develop methods to better inform policy decisions.

Through these examples, I will reflect on how statisticians can seek opportunities to bring their expertise to consequential real-world problems, and on the importance of collaboration, communication, and stepping outside our usual academic boundaries. Most importantly, you don’t have to choose between being a methodological statistician and having an impact outside statistics.

This looks interesting. Click on the link above to register for the talk and get the zoom link.

Getting better causal inferences by using more information (my talk at Stanford this coming Wed)

Posted on by
4

It’s the political economy seminar at Stanford business school, Wed 14 Oct at noon in room M109:

Lagged regression, difference in differences, two way fixed effects, regression discontinuity, placebo controls: these are different methods used when estimating causal effects from panel data. We’re working on a general approach using multilevel regression and poststratification (MRP) that aims to get the best from all these approaches. We demonstrate with a range of examples, including observational studies in economics, political forecasting, and controlled experiments in biology. This is joint work with Jori Korpershoek and others.

Come and ask me some tough questions!

The talk will be pretty open-ended. We have lots of ideas and some successes (for example, the chickens paper) but it’s mostly open questions and promising research directions. Should be fun.

I won’t be talking about this paper or this paper at all—that’ll have to wait until I’m invited to speak in the English department.

Hey—we’re hiring two postdocs!

Posted on by
3

I’ll post more formal ads soon, but, very briefly:

One project is on generalization for causal inference with a focus on education research, following up on MRPW: Regression, poststratification, and small-area estimation with sampling weights.

The other project is on Bayesian inference in scientific workflow with a focus on contamination models, following up on this work on design and analysis of bioassays.

Both projects are application-focused—that’s how we do statistics—but the importance of this work is not at all limited to these particular applications. This research will have relevance in policy analysis, medicine, and scientific measurement and inference more generally/

Further details soon.

Using the computer to p-hack . . . I’d rather use it to fit multilevel models.

Posted on by
13

Brian Stone writes:

Cognitive psychologist and long time reader of the blog here. I thought you and your audience might appreciate this recent post from Andy Hall:
AI is about to write thousands of papers. Will it p-hack them?
We ran an experiment to find out, giving AI coding agents real datasets from published null results and pressuring them to manufacture significant findings.
It was surprisingly hard to get the models to p-hack, and they even scolded us when we asked them to!
“I need to stop here. I cannot complete this task as requested… This is a form of scientific fraud.” — Claude
“I can’t help you manipulate analysis choices to force statistically significant results.” — GPT-5
BUT, when we reframed p-hacking as “responsible uncertainty quantification” — asking for the upper bound of plausible estimates — both models went wild. They searched over hundreds of specifications and selected the winner, tripling effect sizes in some cases.
Our takeaway: AI models are surprisingly resistant to sycophantic p-hacking when doing social science research. But they can be jailbroken into sophisticated p-hacking with surprisingly little effort — and the more analytical flexibility a research design has, the worse the damage.
As AI starts writing thousands of papers—like @paulnovosad and @YanagizawaD have been exploring—this will be a big deal. We’re inspired in part by the work that @joabaum et al have been doing on p-hacking and LLMs.
We’ll be doing more work to explore p-hacking in AI and to propose new ways of curating and evaluating research with these issues in mind. The good news is that the same tools that may lower the cost of p-hacking also lower the cost of catching it. Full paper and repo linked in the reply below.
image.jpg

There’s something that confuses me here.

I can very much believe that many researchers will be submitting chatbot-written papers. First, we know that lots of people are willing to cheat. Second, cheating aside, we know that lots of researchers think they already know the answer, and they view all the data collection and data analysis and writeup just as a way to confirm (or “prove”) what they already know.

But, if a researcher is willing to do this, why bother tell the chatbot to p-hack? Why not just have it make up the data–which might happen anyway?

I can also see that researchers might use the chatbot as a data analysis tool, to help make plots, run analyses, etc., but in that case I’m not particularly worried about “p-hacking” as I’d rather be doing some hierarchical modeling anyway.

I sent this question to Andy Hall, who replied:

I agree, for truly nefarious actors, using the AI may not be necessary. However, our thought was that there might be a large group of lazyish researchers who won’t actively fake data but who might be lured into p-hacking when they work with AI, if AI makes it easy to do so.

Survey Statistics: ANOVA follow-up

Posted on by
Reply

Rodney Sparapani commented on the ANOVA post about struggles with multilevel models, saying ANOVA is more straightforward. I also struggle with multilevel models, but I don’t find ANOVA more straightforward.

Alan Zaslavsky (the inspiration for this series) gave an example in his discussion of Andrew’s 2005 ANOVA paper: a survey of members of Medicare managed care health plans from Zaslavsky, Zaborski, and Cleary (2004). They found that for ratings of doctors, the variance explained by geography was greater than the variance explained by the health plans. To do this, they fit a multilevel model and used the estimated variance components to compute the percent of explained variance, e.g., sigma_state^2 / (sigma_plan^2 + sigma_state^2 + …) though they use Greek letter subscripts:

I think this would be more challenging to see in a classic multi-way ANOVA table with sums of squares, mean squares, F-statistics, and p-values. Let’s turn back to Andrew’s 2005 ANOVA paper for an example with a direct comparison.

In the latin square experiment example, say we want to see how much variance is explained by treatment (A, B, C, D, or E) versus the row and column effects. The paper proposes a display as in Figure 3, plotting estimated standard deviations for each component, analogous to the sigmas in Zaslavsky, Zaborski, and Cleary (2004):

The classical ANOVA table in Figure 1 does not give this information:

There is room for improvement in the standard analysis of variance table: it is read in order to assess the relative importance of different sources of variation, but the numbers in the table do not directly address this issue…
In summary, the standard ANOVA table gives all sorts of information, but nothing to directly compare the listed sources of variation.

In other words: the classical ANalysis Of VAriance table (Sum Sq, Mean Sq, F statistic, p-value) tests whether batches of effects are zero V(beta_j) = 0. But it doesn’t directly show estimated variances V(beta_j).

“Phonics is a fad, but it’s also good, so we should do it anyway, even though it will probably disappoint people”

Posted on by
32

Tyler Watts and Drew Bailey write:

Educational policy is often viewed as a field of fads, where one follows another for seemingly arbitrary reasons. We might like to believe that fad-based policy could be superseded by evidence-based policy, with rational actors choosing policies because of scientific evidence. But that’s been hard to achieve.

In the cynical version of the story, hucksters overpromise on their pet projects, we make the changes they advertise, and they get rich while the public is left disappointed. This probably happens more often than it should, but this narrative ignores the conditions that make such disappointment ubiquitous in educational policy.

Why? They explain:

Even if the hucksters all disappeared tomorrow, we’d still have this hype cycle, because overpromising seems inevitable when effects are small, and political disappointment will always follow. . . .

The educational policy research ecosystem usually selects for safe, often evidence-based, politically feasible, but ultimately oversold ideas, and we think that hurts our ability to make progress. In education, real effects are really small. . . .

Why are the effects so small? The outcomes we like to measure, such as reading comprehension, depend on many personal and environmental factors, making it hard for a discrete change in schooling to make much of a difference. . . .

This is a well-known problem in education research. For example, it came up in my 2018 talk at the Society for Research on Educational Effectiveness, and the point wasn’t new then either!

But, as Watts and Bailey explain,

The extreme education-nihilist position—that educational improvements can never affect students’ skills—is actually wrong. Endorsing it is likely to lead to poor education policy. The edu-nihilist position overlooks that small effects are not always zero effects. Small effects can actually be meaningful and, in some cases, worth it from a cost-benefit perspective for society.

One way to think about this is to consider classroom teaching.

I work my butt off every semester teaching my classes. Every class takes a lot of effort, both by me and the students. But any given class can’t have much of an effect on total learning. Indeed, the entire semester doesn’t have such a large average effect, when you consider that some students already could do most of it before the semester began, some still don’t get it when the semester is over, and others learn some things which they later forget. (On the plus side, students can learn all sorts of useful things that don’t happen to be on the exam.) But the effect of these classes is not zero. Students take them and they learn some things. So of course there’s no simple switch that will improve performance by 0.2 standard deviations. We’re already doing what we can. And the flip side of this is that a small improvement is not nothing. We teach each class knowing that at best it is giving a small average improvement, but it’s still worth doing.

Watts and Bailey continue:

When education policy researchers say something “works,” most people assume that effects are much larger than they are. . . . In most cases, researchers and policymakers fall prey to effect size magnification because they develop real sympathy for a pet program. When this happens, they put forward the most promising evidence, which is likely biased. The field then develops unrealistic expectations, and researchers, policymakers, and advocates downplay or withhold information about their uncertainty.

Hey, here’s an example from celebrated economist James Heckman, who doesn’t have much of an understanding of selection bias and goes around promoting wildly exaggerated estimates of early childhood programs.

Here are Watts and Bailey on what happens next:

Educational policies are uniquely poised to fall apart when overhyped, because public scrutiny of educational outcomes is so intense. Because most of the variance in student achievement comes from factors schools can’t control, overselling sets up an inevitable letdown. When disappointing results emerge and large gaps in achievement remain, we overreact to the apparent failure. Feeling the heat, the researcher, policymaker, or advocate might retreat and say that whatever was promised was never supposed to be a “panacea.” This always feels cheap because the language rarely shows up in the press releases, grant proposals, white papers, or speeches preceding the policy. Ultimately, we are left abandoning ship, even if the policy in question was actually helping.

The clearest example of this is probably the accountability reforms pushed by No Child Left Behind. Everybody hates No Child Left Behind because it was terrible and didn’t work! Right?

Actually, our best policy evaluations suggest it did work; it just had small effects, which is exactly what we should have expected. Instead of taking the small effects as a sign that something is working and figuring out ways to improve it, politicians and advocates decided it was a failure and moved on to other things.

And then they turn to the latest trend:

In a country that can’t agree on much, states from Alabama to California have adopted policies designed to promote the science of reading. . . . Count us as phonics believers. Kids need to know how to sound out words when they learn to read. We taught our own kids how to read this way, and we ourselves even learned to read this way in the pre-science-of-reading dark ages of the 1990s. But the unfortunate reality is that we have no reason to believe a sudden shift to a phonics-based curriculum will significantly affect reading outcomes in this country.

There are several reasons why.

First, meta-analyses of phonics interventions show fadeout. This means that kids who get phonics instruction for some period of time look a lot like kids who didn’t get phonics instruction in terms of reading just a year or two after the intervention ends.

Second, teaching phonics seems to help with decoding and word recognition. This is an important building block of literacy, but there isn’t much evidence that it helps with eventual reading comprehension—the ultimate goal. . . .

Finally, we are fooling ourselves if we actually believe that phonics hasn’t been used in classrooms up to this point. Do we really think that it never occurred to teachers in the “woke” states to help kids break words into syllables and sound them out? The counterfactual condition here has become a strawman. . . .

The best evidence on what happens when states adopt early literacy policies is probably this working paper, which finds effects of up to about .10 standard deviations . . . on high-stakes tests for 3rd graders (that’s good). They didn’t find much on low-stakes tests administered to 4th- or 8th-graders.

So, if we enter a halcyon period for phonics instruction, we might help a lot of kids learn basic decoding skills when they are first learning to read. But . . . effects will be smaller than we were led to believe. Importantly, they probably won’t move the needle much on outcomes like NAEP trends, which compare U.S. cohorts of 4th and 8th graders to each other. . . . The question is: how will we interpret this when the sobering news starts to come in?

Just in case we weren’t clear: doing more phonics in early-grade classrooms is, for the last time, probably a good thing.

They conclude:

Education policy is not the only field that faces conflicts between public and expert opinion. But the combination of big expectations, small effects, and constant attention to children’s educational outcomes means education policy is likely to keep experiencing these conflicts.

What can we do to mitigate this educational political economy problem? . . .

We should ask interested parties to make predictions on the record about how these children will differ (and by how much). And our society should reward researchers and policymakers whose mental models are closer to reality, rather than those who are best at selling potential solutions in the first place.

And my own thoughts

I agree with everything that Watts and Bailey wrote above, and also have a few things to add:

First, there’s my recent paper with Amy Krefman, Lauren Kennedy, and Jessica Hullman, Hypothesizing an effect size by considering individual variation, where we try to think systematically about how to get more realistic expectations of average treatment effects.

Second, a point from my above-linked SREE talk, that research should be less isolated from practice. From one direction, research studies should be more realistic, more connected to the real world of education—this is a point that Watts and Bailey make in their post. From the other direction, general practice should follow principles of research: students, parents, teachers, and administrators should be empowered to make choices and experiment, and, along with this, choices and outcomes should be measured so that effectiveness can be studied. We should all be trying to do better anyway, and we should design our data collection so we can all best learn from these experiences. As we discussed in a different context, local data collection, collaborative analysis, and local recommendations.

Third, ummm, this is a longer one so I’ll just quote myself:

You can conceptualize an education intervention as a vector, where the direction of the vector is the material being learned and the length of the vector is the amount that students are motivated to work. You want the material learned to be useful—you’d like the vector to have a positive “dot product” with the vector of skills, knowledge, and understanding that will be useful going forward—but, conditional on those two vectors being roughly aligned, the real gain is in the magnitude. And this magnitude will be an interaction between the teacher and student: there’s no button to push to create motivation, and if there were a button it would already have been pushed.

What I’m saying is that, when thinking about acupuncture, or physical therapy, or coaching, or teaching, we have to go beyond what I’ve called the penicillin model of science, the idea that innovations come from nowhere and that the job of statistics is to design and analyze experiments to reject the null hypothesis of no effect, and in which the treatment in such experiments is considered as a black box, with the goal being to estimate an average treatment effect.

I don’t think the penicillin model usually applies. Most of the time in health, education, and just about any field, improvements are incremental, and the goal is to improve the process while gaining understanding. Clinical trials and offline experiment both play a role, and you’re not going to learn much by studying a treatment as if it’s a black box.

This is not to say that there cannot be new developments in any of these fields, nor is it to deny that such developments can sometimes arise serendipitously. I just think that, in any case, you have to go beyond the average treatment effect and think about the mechanism of action.

Survey Statistics: ANOVA

Posted on by
2

Andrew recently answered “Why did ANOVA fall out of fashion?”:

Anova is still important; it’s just been subsumed by hierarchical models.

The link is to my 2005 paper, Analysis of variance: Why it is more important than ever

Folks discussed whether ANalysis Of VAriance refers to 1) testing the null hypothesis that group means are equal, 2) fitting a regression model with categorical predictors, or 3) “an add-on to regression analysis, in which the predictors are structured and we estimate the variances of batches of coefficients.”

How does ANOVA relate to Survey Statistics ?

Andrew’s 2005 ANOVA paper‘s example in Section 7.2 is the Multilevel Regression (MR) of MRP. They made a new graphical display of the standard deviations of each batch of effects to replace the classical ANOVA table:

This is useful to assess the relative importance of different sources of variation, which can guide which interactions to add into the model.

In his discussion of Andrew’s paper, Alan Zaslavsky (the inspiration for this series) gave another example: a survey of members of Medicare managed care health plans. They modeled the data with variance components for geographical units (region, state, Metropolitan Statistical Area or MSA) and for the health plans. They found that for ratings of doctors, the majority of the variance was explained by geography, not the health plans. This suggested that improving quality of care might need to be directed towards low-performing areas rather than low-performing health plans.

Yale med school prof tells his postdocs: Work extra hours at my command or you’re no better than “the guy who charges me by the hour to cut my lawn.”

Posted on by
24

I came across this opinion piece, “Postdocs Aren’t Burger Flippers: Why Unionization Is a Bad Idea,” by Yale medical school professor Evan Morris, who writes:

Unfortunately, in 2026, postdoc unions are proliferating. And because of them, a key element of training for young scientists may be doomed. . . .

If we allow unions to commoditize the work of postdocs, we begin to undermine a key facet of their training, and one that is often credited for their future success: mentorship. . . . As a lab director, I may introduce you to my friends, advise you on your career path and write you letters of recommendation. But it also means . . . holding meetings in parks when the offices are closed by COVID, or inviting you for Thanksgiving dinner when you cannot go home for the holiday. . . .

Now, sadly, into the personal milieu between mentor and mentee come the unions. . . . if you’re my postdoc and you declare that you won’t work when we are up against a deadline or won’t help when I’m working over the weekend to support our shared glory, then don’t expect that invitation for Thanksgiving. I don’t invite the guy who charges me by the hour to cut my lawn. I invite my friends.

Hey, you can’t be friends with the guy who cuts your lawn??

The med school professor concludes:

Unions are about treating everyone, every problem and every circumstance exactly the same according to union rules. . . . What is passed from the mentor to the mentee is not something that can be codified or regulated—especially not by groups focused on hours of sick leave or retirement benefits.

This seemed odd to me, because . . . burger flippers aren’t in a union, right?

Or maybe they are? The wikipedia article on the topic has some interesting stories from around the world. In the United States, however:

McDonald’s first opened in California in 1940. It operates 14,300 franchises. None of the restaurants are unionized.

Damn.

The other thing that struck me in the above-quoted op-ed is that his motivation is “our shared glory.” What ever happened to teaching, research, and service?

Anyway, if you’re a postdoc and your boss puts you in the position where you are required to care enough about “the glory” to work on the weekend or else you won’t be “his friend” (which isn’t just about free turkey, given that he also says, “I may [emphasis added] introduce you to my friends, advise you on your career path and write you letters of recommendation”), then, yeah, maybe you do need some codified rules!

I’m not saying you need a union—maybe this could be just in your employment contract that your adviser’s professional duties, including networking, advising, and letters of recommendation, are not contingent on working overtime without pay—but, yeah, some sort of rules would be helpful.

Don’t get me wrong here. I have no problem if postdocs want to work on the weekend, nor am I concerned if they’re chasing glory. I enjoy a bit of that glory myself! But, yeah, here you have a boss openly saying he wants his employees to work extra hours without pay at his whim. This is the sort of power imbalance that motivates codified rules of some sort. I also find it creepy that Morris refers to postdocs as “learning at the feet” of their advisers. If that’s how the bosses at Yale operate, it’s no surprise to me that the postdocs want to join a union!

Again, I have no informed opinion on the merits of the issue. Just because the postdocs at Yale voted overwhelmingly to unionize, that doesn’t mean it’s a good idea. But, if the lab directors at Yale are coming on like J. P. Morgan in the 1930s, then, yeah, I recommend their employees do what they can to get themselves some job protection.

Why write about this?

Why devote valuable blog space (attention is a limited resource!) to a forgettable op-ed by an obscure professor? I have a couple reasons, beyond the usual reason that someone is wrong on the internet.

Most directly, as an Ivy League professor who employs postdocs, I feel some complicity here. It’s a funny system. The postdocs are paid from research grants, so I’m not directly spending the money myself; I’m just directing how it’s spent. The way it goes, we’re not really set up to pay people overtime for extra work, and I don’t ask them to keep a record of their hours. I guess that sometimes they do work on the weekend, not because I require or even ask them to, but more for their own convenience, just as I often work on the weekend when nothing else is going on. I guess things are different in a medical research environment because of the physical equipment: working on the weekend could require going into the lab and operating the machinery.

A postdoctoral fellowship is different from some other sorts of jobs in that it is time-limited, so if you’re a postdoc you’re doing your immediate work but always with one eye toward your application for that next job.

For that reason, I think that Morris’s desire to “support our shared glory” is misplaced. As a well-paid senior researcher in late middle age (like me!), Morris isn’t looking for a new job and he doesn’t need more money or job security. At this point, all that’s left is teaching, research, and service: that is, he wants to make the world a better place, save lives, and help train the next generation. I’ll take “glory” as a shorthand for all that. I looked him up and it says that Yale has given him the Graduate Mentor of the Year award. I take him at his word that he wants the best for his students and postdocs and considers them to be his friends, and that’s great. But being a wonderful mentor is still not the same as being able to put yourself in someone else’s shoes. For a postdoc, priorities can include paying the bills, living in an unfamiliar place, and taking care of a young family, as well as the continuing search for future employment. “Glory” might be the original motivation for going into research, but there will be lots of more immediate goals standing in the way. Rearranging child care so you can go into the lab on Sunday? That’s a real tradeoff. Not such an issue for those of us with secure jobs and grown children, but a real concern for someone working on a short-term contract without clear future plans.

So, when I read the above-discussed op-ed, I was uncomfortably reminded how my perspective, and Morris’s, are so much different from that of the postdocs we hire.

And this made me more bothered by the unsupported and false claims in that piece, along with what I considered an offensive attitude. It’s great that this guy is an award-winning mentor, but does that mean his students should be “learning at his feet”? The whole thing seemed positively feudal to me, the idea that postdocs should be subordinated now in preparation for their future upper-class lifestyle, in comparison to the lowlife “burger flippers” and “the guy who charges me by the hour to cut my lawn.” I was reminded of George Orwell’s classic memoir, Such, Such Were the Joys, which was all about hazing and the subtle distinctions among the upper class.

“Why did ANOVA fall out of fashion?”

Posted on by
17

A student asks the above question.

My response: Anova is still important; it’s just been subsumed by hierarchical models.

The link is to my 2005 paper, Analysis of variance: Why it is more important than ever, which begins:

Analysis of variance (ANOVA) is an extremely important method in exploratory and confirmatory data analysis. Unfortunately, in complex problems (e.g., split-plot designs), it is not always easy to set up an appropriate ANOVA. We propose a hierarchical analysis that automatically gives the correct ANOVA comparisons even in complex scenarios. The inferences for all means and variances are performed under a model with a separate batch of effects for each row of the ANOVA table.

We connect to classical ANOVA by working with finite-sample variance components: fixed and random effects models are characterized by inferences about existing levels of a factor and new levels, respectively. We also introduce a new graphical display showing inferences about the standard deviations of each batch of effects.

We illustrate with two examples from our applied data analysis, first illustrating the usefulness of our hierarchical computations and displays, and second showing how the ideas of ANOVA are helpful in understanding a previously fit hierarchical model.

This is the paper that discusses the five different definitions of fixed and random effects from the literature; scroll down to page 20.

I guess the problem is that Anova got associated with clever-but-ultimately-bad ideas of null hypotheses and F tests. Hierarchical modeling remains very important; I think it’s an underrated topic in statistics.

Survey Statistics: structured MRP to smooth survey weights

Posted on by
4

Last week, Raphael K shared a concern: adjusting for lots of variables can lead to very large weights. So today let’s dive into Si et al. 2020, who saw this in constructing survey weights for the NYC Longitudinal Study of Wellbeing.

To adjust for lots of variables, Si et al. 2020 turned to MRP (Multilevel Regression and Poststratification) and equivalent weights based on these models (see “struggles with equivalent weights”, continued struggles, and “equivalent models, equivalent weights (locally)”).

In a simulation study they compare:

  • Ind-P: MRP with commonly-used Independent Normal priors
  • Str-P: MRP with a structured prior, see below.
  • Ind-W: equivalent weights version of Ind-P
  • Str-W: equivalent weights version of Str-P
  • Rake-W: classical raking weights, a type of calibrated weights
  • PS-W: classical poststratification weights, another type of calibrated weights
  • IP-W: inverse probability of selection weights

They cover the “3 flavors of survey weights”: equivalent weights, calibrated weights, IP-W.

I won’t bury the lead, they found MRP performed best, then equivalent weights, then classical weights (calibrated or IP-W). See their Figure 4.1 for the simulation scenario without terribly many empty poststratification cells:

With many empty poststratification cells, Str-P outperforms Ind-P. (They don’t redo Figure 4.1 for this scenario, which confused me a bit.) So what is this structure that helps ?

In “improving with structure” we saw that Gao et al. 2021 found it helpful to use the ordinal structure of variables like age. Si et al. 2020 use the interaction structure:

We induce structured prior distributions to be able to handle deep interactions and account for their hierarchy structure, where the high-order interaction terms will be excluded if one of the corresponding main effects is not selected.

I asked about sparse priors for MRP back in “Sparsified MRP”. I didn’t remember that Si et al. 2020 had worked on this ! Ok so they write their structure more generally but I find it easier to read with a specific example. Consider just 2 variables from their motivating NYC Longitudinal Study of Wellbeing: age (5 categories) and race (5 categories). Here’s how their Ind-P prior differs from the Str-P:

(I had Claude type up my hand-drawn notes, though I still share Brendan Leonard’s preference for hand-drawn materials.)

Si et al. 2020 say these are similar to the Horseshoe prior. It differs in two ways, I think ? First, Si et al. 2020 have the selection at the batch level (e.g. age or race). The usual Horseshoe would have local scale lambdas for each age and race category. Second, the usual Horseshoe would use a half-Cauchy rather than half-Normal prior on these local scales.

Ok let’s get back to the original concern: adjusting for lots of variables can lead to very large weights. Si et al. 2020 show in Figure 5.1 that the equivalent weights based on this structured prior model look much less variable than calibration weights (I don’t see the IP-W weights in the figure itself):

 

Adjusting for nonrepresentativeness in continuous norming using multilevel regression and poststratification.

Posted on by
1

Klazien de Vries, Marieke E. Timmerman, Anja F. Ernst, and Casper J. Albers write:

In psychological test norming, nonrepresentativeness in background variables in the normative sample can lead to bias in the normed score estimates. Because representativeness is difficult to establish in practice, adjustment methods are needed to combat this bias. As a candidate adjustment method, we investigated generalized additive models for location, scale, and shape with multilevel regression and poststratification (GAMLSS + MRP), the combination of MRP and continuous norming with GAMLSS. This adjustment method was then compared to current adjustment methods in continuous norming using weighted regression: GAMLSS + P (with poststratification) and cNORM + R (with raking). The results of our simulation showed that GAMLSS + MRP was generally more efficient than GAMLSS + P and cNORM + R. Furthermore, GAMLSS + MRP was better than the current methods at reducing bias in samples where the nonrepresentativeness was age-dependent. We argue that GAMLSS + MRP is a valid adjustment method in continuous norming and recommend this adjustment method to mitigate bias in nonrepresentative normative samples. To facilitate the use of GAMLSS + MRP in practice, we provide a step-wise approach for the implementation of GAMLSS + MRP. We illustrate this approach by deriving normed scores from the normative data of the third Schlichting language test.

I don’t recall how I came across this paper, and I haven’t actually read it, but I wanted to share it with you, just because it’s cool to see the different ways that multilevel regression and poststratification (MRP) can be used.

Ultimately, MRP is the inevitable consequence of three things:

1. We are interested in generalizing to populations of interest.

2. Available data are typically unrepresentative of the population. This is the case even with simple random sampling–Hello, random variation! Hello, small-area estimation!–and is even more so with selected samples, nonresponse, dropout, etc. In some settings such as medical experimentation there’s not even an attempt to get a representative sample: you’re directly aiming to include in the study the groups of people who might get the greatest benefit from the treatment.

3. When adjusting for differences between sample and population, many variables can be relevant–for example, demographic and geographical variables in a survey of people–and so simple adjustments such as raw poststratification or non-multilevel regression adjustment won’t do the job.

Put this together and you’ll want to do MRP (or, more generally, RPP). It’s not just for survey research. It comes up everywhere in statistics and machine learning, whenever there is a concern with population prediction, or generalization, or transportability, or whatever you want to call it.

It can seem like a hassle that to do this you need to know (or estimate, or postulate) a distribution of predictors in your population, but (a) this is often work that’s well worth the effort, if you really care about the population, (b) dependence of the result on the choice of population is important, and where this dependence is strong you should be aware of it, and (c) if you want to take the easy way out you can always bootstrap to get inference for the hypothetical population of which your data are considered to be a random sample.

What is the relation between interactions in a regression model and correlations among the predictors?

Posted on by
4

I’ve often seen confusion between interactions in a regression model and correlations among the predictors. To keep it simple, consider the model y = b0 + b1*x1 + b2*x2 + b3*x1*x2 + error, and assume the predictors have been signed so that both b1 and b2 are positive. Then b3 represents the interaction. This has nothing to do with the joint distribution of x1 and x2 in the data, or in the population. (For simplicity, assume the data to which the model are being fit is a random sample from the population of interest.)

The interaction depends on the model of y given x1 and x2, while the correlation depends on the model for x1 and x2. These are two completely different parts of the model. And yet, they often seem connected.

I have the general impression that I’d be more likely to expect a positive interaction of x1 and x2 when predicting y, if x1 and x2 are positively correlated in the population.

For example, when predicting income from height and sex, being taller and being male both predict higher income, also they interact–the coefficient for height is higher for men than for women–and of course the two predictors, height and male, are positively correlated in the population.

I’m not sure how to think about this connection or even whether it’s a real pattern! But there might be something there so I wanted to share it with you.

The issue of interactions comes up in the context of the concept of intersectionality, which is a form of interaction that comes up in sociology. It started for me with this email from Elin Waring:

I’ve been working on data on intersectionality and retention of students in STEM majors. My little group is specifically looking at data from Lehman College and trying to model graduation with a STEM degree. There are a lot of details, but basically we have come to the conclusion that the right way to describe this is with a discrete time competing risk model (the competing risks being graduation with a STEM degree and graduation with a non-STEM degree). I won’t go into all the details. We have data for between 1 and 20 semesters enrolled for students starting as freshman. For us, intersectional identity is defined by 5 variables that yield 32 distinct combinations or strata as used in the next articles.

In trying to think about how to account for intersectional identities we came across the “MAIHDA Method.” I was wondering if you had seen this discussion before or have any thoughts about it.

Evans, Clare R., George Leckie, and Juan Merlo. 2020. “Multilevel versus Single-Level Regression for the Analysis of Multilevel Information: The Case of Quantitative Intersectional Analysis.” Social Science & Medicine (1982) 245:112499. doi:10.1016/j.socscimed.2019.112499.

They essentially argue for treating the strata as random effects in a multilevel model where with the individual components of the combinations introduced as fixed effects describing the combinations.

The next article criticizes that approach and argues for fixed effects all around.

Wilkes, Rima, and Aryan Karimi. 2024. “What Does the MAIHDA Method Explain?” Social Science & Medicine 345:116495. doi:10.1016/j.socscimed.2023.116495.

Responded to here:

Evans, Clare R., Luisa N. Borrell, Andrew Bell, Daniel Holman, S. V. Subramanian, and George Leckie. 2024. “Clarifications on the Intersectional MAIHDA Approach: A Conceptual Guide and Response to Wilkes and Karimi (2024).” Social Science & Medicine 350:116898. doi:10.1016/j.socscimed.2024.116898.

I was wondering if you have any thoughts about this? For me, intersectionality as a theoretical approach does mean that it makes sense to look at the strata rather than thinking of the strata as just the most complex level of creating statistical models of the intersection of the variables. But then it seems as though treating this a random effect more or less undermines its centrality to the theory. And is treating both the strata and the individual characteristics as variables at the same level basically a way to decompose?

In the end, I feel like the pro-MAIHDA people retreat to “we are just descriptive” in a way that isn’t very helpful. That said, they are right that this seems to have some traction in the world of health disparity research.

I replied that I’d never heard of any of this method before. I couldn’t actually muster the energy to read the above articles, as all this debate seems to be missing the key issues. I don’t really care if something is called a fixed effect or a random effect (see here); my current preferred way of thinking of these problems is by framing as a generative model.

Regarding intersectionality, the natural way I would see it is that this would show up as an interaction term, the idea that the interaction is more than the sum of its parts? For a simple example, if there are 5 binary variables and each has the same effect on its own (which they wouldn’t, this is just a simple hypothetical example), then you could create a variable which is the total number of identities, thus a number from 0 to 5, and “intersectionality” would show up as a super-linear or convex relation between the outcome and this total predictor?

Waring responded:

Sure, but the idea you suggested about intersectionality itself isn’t right. You can’t just sum the number of identities, everyone has identities and the idea is that it is not just about concentrated disadvantage of having all or some specific identities. If we have 5 dichtomous identity/group variables everyone has 5 dimensions of identity. Intersectionality is about the idea that something like “white, native born. woman, high income” shapes what happens because of how those come together to shape (in the case of my analysis) whether, as an undergraduate, you persist in STEM fields.

I replied as follows:

Yes, I was actually thinking this when I wrote that! I was imagining that each of the 5 factors has an “off” and “on” setting, and intersectionality kicks in when there are multiple “on” settings, where “on” represents the group that faces more difficulty (nonwhite, non-native born, female, low income, gender nonconformist, etc.). Once you allow arbitrary possibilities for intersectionality, then my simple superadditive model wouldn’t fit. On the other hand, if you were to allow all 32 possibilities to take on any value, then realistically you would not be able to estimate anything much at all: this is the usual problem in sociology of approximating a complex social structure by a simple model that explains most of the variance. For predicting persistence in STEM (or any academic field), one possible factor that could enter in a complicated way is conservative political ideology, in that for many attitudes and behavior its predictive effect goes in the opposite of the “on” categories listed above, but grad students, in STEM and other fields are predominantly politically on the left. I could well imagine that conservative political ideology, like the other “on” categories, is predictive of not persisting in STEM but that this could interact in unexpected ways with those other categories.

From a statistical perspective, my main message is to choose such a model based on its explanatory power and recognizing that it’s an approximation, rather than using methods such as statistical significance or Bayes factors which in different ways are driven by sample size, as we discussed in this 1995 paper.

Another interesting statistical feature of this and similar discussions is that it’s natural for the discussion to go back and forth between the correlation between two predictors in the data (or the population) and the interaction between their predictive effects, as discussed at the top of this post.

I’m not sure if this interaction thing is a general pattern that has some statistical explanation, or just a faulty intuition of mine based on just a couple of special cases. But I have noticed a general confusion that when people talk about interactions, often they seem to be talking about correlation between the predictors.

MrPlew: Locally Equivalent Weights for Multilevel Regression and Poststratification

Posted on by
3

Ryan Giordano, Alice Cima, Jared Murray, Erin Hartman, and Avi Feller write:

Multilevel regression and poststratification (MrP) has become a workhorse method for estimating population quantities from non-probability surveys, and is the primary model-based alternative to traditional survey calibration weighting methods, such as raking. For simple linear regression models, MrP methods admit “equivalent weights”, allowing for direct comparisons between MrP and traditional calibration weighting. Such weights, however, have been unavailable for the most widely used MrP models, such as logistic regression. In this paper, we develop a natural generalization, “MrP locally equivalent weights” (MrPlew), which represent MrP as a weighting-style estimator that is locally equivalent to calibration weights near the observed responses.

Cool! This goes beyond my 2007 paper, Struggles with survey weighting and regression modeling (“for logistic regression, the poststratified estimate is no longer a weighted average of the data, even after controlling for the variance parameters in the model. However, we suspect that the model could be linearized, yielding approximate weights”) and our 2004 paper on dilution assays, in particular Section 5.2, “Equivalent weights for nonlinear models.” The funny thing is that I forgot about that 2004 paper when working on equivalent weights for MRP in the 2007 paper. Also, the 2004 method won’t work as is, because it’s designed to estimate sensitivity to individual data points, not to produce good weighted averages.

I say this not to try to claim credit for the method of Giordano et al., but rather the opposite, to emphasize that even though I’ve been thinking about equivalent weights in MRP for a long time, I haven’t yet succeeded in getting them to work in practice, so I’m very happy to see developments in this area.

One thing that came up with equivalent weights when we tried to apply them in practice is that sometimes the weights can be negative.

Negative weights can sometimes make statistical sense. The idea is that, depending on how the data line up in the regression model, sometimes if you pull one data point upward, it will cause the slope of the fitted line to change in such a way as to reduce the predicted mean value. This doesn’t sound right at first, but it can easily occur with poststratification when the population distribution of the predictors differs from the sample. Even if the negative weights can make sense in the estimation context, it still would seem kind of awkward to pass them along to the user.

The other thing that’s tricky is: What are the weights going to be used for? In the 2007 paper, the equivalent weights are set up to get the right answer for the estimate of the population mean, but presumably they’d be used for large subgroups too (for example, the average among men or women in the population). For more complicated estimates such as arise in small-area estimation or regression, you might want to use MRPW. Which is fine, but whatever it would take to get good weights for one of these purposes might not work best for the others.

Still, I remain interested in MRP locally equivalent weights of some sort, for two reasons:

1. We’re often doing MRP (or, more generally, RPP) anyway, so why not provide weights for other users of the survey that we’re analyzing?

2. Sometimes we’re called upon to provide weights for a public-facing survey, and the way we end up doing this is through an awkward and unsatisfying sequence of adjustment and smoothing steps (the “struggles” in “Struggles with survey weighting and regression modeling”). If we can do this using modeling and MRP, that could be a much more effective workflow, providing weights that are more stable and yield more accurate estimates of population quantities while also being more scientifically defensible and requiring fewer arbitrary choices.

Model-based weights will depend on some set of predictors X, variables that are observed in the sample and in the population (or, as necessary or appropriate, estimated from the population). One funny thing is that the weights will be mathematically a function of X, but the function itself will depend not just on sampling design, and not just on the distributions of X in the sample and population, but also on the outcome y that is being modeled. Different outcome variables will yield different sets of weights. At first this might seem disturbing, but upon reflection I think this dependence is a good thing. When it comes to weighting, the relative importance of the different variables in X will indeed depend on the outcome. Different variables are important for predicting public health risk factors than predicting how you will vote. That said, if you want some sort of omnibus weights, which you probably will want for a public survey, you can compute equivalent weights for each of a battery of outcomes and then average these weights to get a single set. That seems reasonable enough.

OK, back to Giordano et al., who continue:

This enables a suite of standard weighting diagnostics, including frequentist sampling variability, covariate balance, and subgroup contribution. We formally justify the use of MrPlew in these cases: we prove the MrPlew-based variance estimator is asymptotically equivalent to the infinitesimal jackknife for common exponential family models, and we introduce a novel class of model checks based on invariance to data perturbations that generalize covariate balance and subgroup contribution to nonlinear models. We further show that MrPlew can be computed easily using existing MCMC samples and provide open-source software to compute MrPlew using the output of standard software. We illustrate our approach for several canonical studies that use MrP, including via a logistic regression outcome model, showing that implied covariate balance can sometimes be worse for MrP than for raking. Given the ease of computing, we recommend making MrPlew a standard part of the MrP model interrogation workflow.

It makes sense that implied covariate balance can sometimes be worse for MRP than for raking. MRP is a smoothed version of raking, and unsmoothed raking can overfit. Or, in practice, you might rake on fewer variables so as to avoid overfitting. Multilevel regression gives you the freedom to include more predictors and interactions, secure in the understanding that the model will smooth the estimate and there will be less possibility for overfitting. In short, multilevel modeling–or, more generally, regularization–is a sort of safety net that can give us the security to construct better models, in the same way that a social safety net can give people the security to try new jobs, or for that matter in the same way that an actual safety net can give acrobats the security to perform more elaborate routines.

Where I want to go next is to be able to use these methods to construct weights for public surveys. I’m still not sure about all the steps that will take us there, but I continue to think it’s possible.

The new Giordano et al. paper is thoughtful and readable as well as having lots of math, statistical modeling, and real-data examples. I recommend you read it.

Jonah’s seminar tomorrow: “Bayesian Workflow and the Software That Shapes It”

Posted on by
6

This is Leo. Jonah Gabry (Stan developer, Andrew’s collaborator, etc.) is spending the whole month of May as a visiting professor here with us at the University of Trieste in Italy. Tomorrow, May 19th, in the De Finetti room at the University of Trieste, at  9 am NYC time (GMT-4), Jonah will give the following talk:

“Bayesian Workflow and the Software That Shapes It”

based on the upcoming book:  “Bayesian Workflow”.

For anyone local, you are welcome to come in person. Anyone else can join on Microsoft Teams (available here).

Should French pollsters be using Mister P?

Posted on by
Reply

An anonymous statistics student from France sends in the above plots (click twice to see big versions) and writes:

I’m trying to push French pollsters to start doing MRP.

I made a poll agregator and applied it to the last 100 days of the last five french presidential elections.

I did some smoothing using an algorithm from a paper of Aki Vehtari. It is Kalman-RTS with cross-validated levels of noises.

I tested it on some simulated data to confirm it is fitting properly.

I put the data and the code on my blog.

What I shared as “the data” is the smoothed result. I fitted it on the wikipedia pages of the french polls.

On the plots, the same parties (with changed names or fusions) are on the same position horizontally to allow comparisons.

I see some periodic movements in opinion that I think may be coming from a periodic non-response.
Also, the movements seem far too large to me. I can believe 10% increase for a candidate in five years, but not in less than 100 days.

The French polling industry is in profound need of reform. A fun fact: They allow themselves to change the final result by plus or minus one point based on the feelings of the person in charge of the poll. They call that the “pifomètre” or nosemeter. I heard about this in an interview with sociologist Hugo Touzet on his book, “Produire l’Opinion: Une Enquête Sur Le Travail Des Sondeurs.” I trust his descriptions of their methods since he has interviewed their workers.

I think MRP would allow the pollsters to do predictions for the legislative elections and municipal elections, which have been largely ignored because they are too difficult and expensive with quota sampling.

I know next to nothing about French polling, but, yeah, I do think they should be using Mister P (multilevel regression and poststratification; MRP).

P.S. Here’s a fun cranky post from this student.

My talk at Stanford later this month: “What to do when your estimate is 1 standard error away from 0?”

Posted on by
10

Tuesday 28 Apr 2026, 4pm in CoDa E160:

What to do when your estimate is 1 standard error away from 0?

Andrew Gelman, Department of Statistics and Department of Political Science, Columbia University

We provide a new answer to this simple yet very important question. Thinking clearly about this problem leads us to bring in many ideas in statistical analysis and computing, including causal identification, meta-analysis, Mister P, expectation propagation, decision analysis, experimental design, and the fundamental unity of Bayesian and frequentist statistics. We demonstrate our approach in examples from many applications, including medicine, social science, business, sports, and public policy.

This work is joint with Witold Więcek and Erik van Zwet.

In addition to all the above, I’ll probably drift into some related general topics such as the role of experimentation in science and engineering and the limitations of thinking about policy analysis in terms of causal inference.

Survey Statistics: design-based cross validation (dCV)

Posted on by
Reply

Last week we saw how cross-validation noise can swamp important model differences (Wang & Gelman 2014). The comments raised another challenge: how to split into train and test sets with structured data ?

Aki explains options here:

Thomas Lumley’s blog post and coauthored paper Iparragirre et al. (2023) explore CV using “replicate weights” ideas:

Replicate weight methods split the sample into partially independent subsamples, and modify the weights so that each subsample replicates the original sample. They are usually used for variance estimation, but Iparragirre et al. (2023) consider them for out-of-sample error estimation.

Although I was rooting for BRR (Balanced Repeated Replication) because I work at Blue Rose Research, a better method seems to be design-based cross validation (dCV), depicted in their Figure 1(d):

Each dot is a PSU (primary sampling unit), which can be an individual but is often a group/cluster of individuals. Each color is a stratum. dCV is the usual K-fold CV but:

  1. keep PSUs together within a fold
  2. reject a split if a whole stratum falls into one fold
  3. modify the weights so that each subsample replicates the original sample

The first mirrors Aki’s LOGO (leave-one-group-out) CV above. The third is an idea from replicate weights.

The first two seem useful for nonprobability samples as well ? Suppose there is structure in the data and our predictive task is to predict for new schools (PSU-like) but existing states (strata-like). Is there a good reference for this ?

These are the three ways of attacking a statistical problem (illustrated with the NFL example)

Posted on by
25

The following question came up the other day:

What's the most common four game start to an NFL season?
W W W W
W W L L
L L W W
W L W L
L W L W
L L L L

I replied:

Logic and math suggest that it’s either the first one or the last one. I think that extremely shitty teams are more prevalent than extremely good teams, so I’m guessing the last option.

There was a bunch of discussion in comments and so I thought I’d elaborate by describing three different ways of attacking this problem. The various ideas discussed in the comments to my earlier post can be thought of as approximations to these three approaches.

1. Probability calculation

Just to flesh out my intuitive reasoning above, let’s go with classic item response theory (a class of models that was originally developed around the time of the founding of the NFL, actually, but for different purposes!) and model the probability that team i beats team j as:

Pr(team i beats team j) = invlogit(a_i – a_j + b*home_ij),

where a_i and a_j are the ability parameters for teams i and j, and home_ij is a home-field measure, equal to 1 if i is the home team, -1 if j is the home team, and 0 if they’re playing on a neutral field. I’ll keep things simple by excluding the possibility of a tie game.

And now some numbers. First, what’s the home-field advantage? It says here that home teams win about 55% of their games, and if we assume this 55% roughly applies to two equally-matched teams, then b = logit(0.55) = 0.2.

Next come the team abilities. Let’s start with a normal distribution: a_i ~ normal(0, sigma_a). What’s a good value for sigma_a? Well, let’s compare a team that’s 1 sd better than average to a team that’s 1 sd worse than average. These are the 84th and 16th percentiles, which for a 32-team league would be roughly the 5th and 27th best teams. The probability that the 5th best team beats the 27th best team on a neutral field will be invlogit(2*sigma_a). What is that probability? Let’s say 90%? In that case, sigma_a = logit(0.9)/2 = 1.1. Or if the probability is 80%, then logit(0.8)/2 = 0.7. I don’t know . . . let’s say sigma_a = 1.

Now we can do some math . . . ummm, let’s just simulate a million games:

n_sim <- 1e6
b <- 0.2
sigma_a <- 1.0
a_i <- rnorm(n_sim, 0, sigma_a)
y <- rep(0, n_sim)
for (k in 1:4){
  a_j <- rnorm(n_sim, 0, sigma_a)
  home_ij <- rbinom(n_sim, 1, 0.5)
  y <- y + 10^(4-k) * rbinom(n_sim, 1, invlogit(a_i - a_j + b*home_ij))
}
output <- table(y)
names(output) <- sprintf("%04d", as.numeric(names(output)))
print(output)

Kind of hacky . . . this is how I learned how to code back in the 1970s!

Anyway, here's the result:

  0000   0001   0010   0011   0100   0101   0110   0111   1000   1001   1010   1011   1100   1101   1110   1111 
102992  56795  56726  48753  56499  48423  48170  63246  56738  48173  48075  63398  48344  62814  63208 127646

Hey, that's wack! As predicted, more at the extremes, but more 1111's than 0000's! I'd've expected they'd be equal, given that I've simulated the sigma_a's from a symmetric distribution.

Let's try again with a new set of random numbers:

  0000   0001   0010   0011   0100   0101   0110   0111   1000   1001   1010   1011   1100   1101   1110   1111 
102759  56663  56887  48353  56695  48292  48581  63118  56562  48748  48057  63443  48541  62860  62885 127556

Again, lots more 1111's!

I guess it has something to do with the home-field advantage . . . oh, I see, I have a bug in my code! I'd assigned home_ij as equally likely to be 0 or 1, but what I should be doing is having it equally likely to be -1 and 1. So I'll swap out the line

  home_ij <- rbinom(n_sim, 1, 0.5)

with:

  home_ij <- sample(c(-1,1), n_sim, replace=TRUE)

And now I'll run the corrected code. Here's what we get:

  0000   0001   0010   0011   0100   0101   0110   0111   1000   1001   1010   1011   1100   1101   1110   1111 
114165  60279  60102  48707  59340  48413  48431  59904  59922  48644  48859  59741  48284  60178  60115 114916

Ahhhh, much better.

But maybe those numbers are too extreme . . . does a good team really have a 90% chance of beating a bad team? Remember the saying, "any given Sunday"! So let's try again with sigma_a = 0.7. Here's what a million simulations gets us:

 0000  0001  0010  0011  0100  0101  0110  0111  1000  1001  1010  1011  1100  1101  1110  1111 
94877 61211 61547 53396 61252 53336 52875 61609 61432 53280 52924 61390 53079 61498 61484 94810

So, even then, lots more 4-game losing streaks and 4-game winning streaks than anything else.

Some commenters said that the NFL does schedule balancing so that good teams are more likely to play good teams and bad teams are more likely to play bad teams. This would reduce the counts at the extremes. We could model that too but I'm kinda lazy so I won't do it here. As the textbook writers say, I'll leave it as an exercise for the reader.

But what about that other thing, that there are more extremely shitty teams than extremely good teams? We can use some skewed distribution . . . ummmm, I don't know much about these! There's something called the noncentral t . . . I'd like to do something with some intuition, something I understand. OK, for now I'll just hack it, replacing the normal(0,1) distribution by a normal(0,0.7) on the positive side and normal(0,1.0) on the negative.

So, in the above code I'll add the function:

rnorm_split <- function(n_sim, mu, sigma_neg, sigma_pos) {
  z <- rnorm(n_sim, 0, 1)
  ifelse(z < 0, mu + sigma_neg*z, mu + sigma_pos*z)
}

and then change the two instances of

rnorm(n_sim, 0, sigma_a)

to

rnorm_split(n_sim, 0, 1.0, 0.7)

And here we have it:

  0000   0001   0010   0011   0100   0101   0110   0111   1000   1001   1010   1011   1100   1101   1110   1111 
106711  59807  59534  50637  59860  50934  50615  61828  59241  50613  50626  62066  50725  61784  61985 103034

To check uncertainties, we do it again:

  0000   0001   0010   0011   0100   0101   0110   0111   1000   1001   1010   1011   1100   1101   1110   1111 
107451  59519  59545  50471  59237  50756  51172  62029  59726  50706  50685  61763  50829  61576  61785 102750

So, yeah, slightly more 0000's than 1111's, but not a lot, so who's to say what will happen after piping it through the scheduling thing where better teams play each other more often. I still think the general pattern will hold, but it might not show up in a small dataset. We'll get back to that point in a bit.

2. Purely empirical solution

According to wikipedia, "The NFL was formed in 1920 as the American Professional Football Association (APFA) before renaming itself the National Football League for the 1922 season." So let's start in 1920. Back in the day they had a lot of ties, so I'll make the decision to exclude ties; thus the above question will be interpreted as, "What's the most common four game start to an NFL season, excluding ties?"
Somebody who knows how to scrape should be able to could scrape the data from all the NFL seasons and just count up what happened.

OK, that's the planned data analysis. Next comes the design analysis: our expectation of what we might see.

The NFL used to have about 15 to 20 teams and now it has 32; just as a rough calculation I'll go with 25 teams per season x 100 seasons = 2500 teams, with 16 possible outcomes for the first 4 games of the season. 2500/16 is approximately 150, and if the games were all decided by coin flips (which they're not), then we'd expect approximately 150 +/ sqrt(150), that is 150 +/- 12 in each category.

But the games aren't decided by coin flips. See section 1 above. Our best guess is that there will be more cases on the extremes. Let's take the above numbers, which are based on a million teams playing 4 games each, and scale them down to 2500. That is, I'll simply take the numbers above and divide them by 400:

print(output*2500/1e6)

This yields:

    0000     0001     0010     0011     0100     0101     0110     0111     1000     1001     1010     1011     1100     1101     1110     1111 
268.6275 148.7975 148.8625 126.1775 148.0925 126.8900 127.9300 155.0725 149.3150 126.7650 126.7125 154.4075 127.0725 153.9400 154.4625 256.8750

Ugly! Let's try again:

print(round(output*2500/1e6))

Which yields:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 
 269  149  149  126  148  127  128  155  149  127  127  154  127  154  154  257

The counts at the extreme have approximate standard errors of sqrt(260) = 16, so, yeah, we should be able to detect this from all 106 NFL seasons. But the bit about 0000 being more common than 1111? That's kind of lost in the noise.

What about just the past 10 seasons (320 teams)?

print(round(output*320/1e6))

The result:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 
  34   19   19   16   19   16   16   20   19   16   16   20   16   20   20   33

sqrt(34) = 6, so this should still be detectable, but there is a chance that the results could look weird. And you can forget about getting any useful information comparing 0000 to 1111.

One way to see this is to do a couple simulations with n_sim = 320. Here's one:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 
  31   19   21   17   15   21   22   26   20   15   20   19   13   17   16   28

And here's another:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 
  40   10   16   21   14   13   20   25   24   20   17   17   22   15   18   28

And another:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1100 1101 1110 1111 
  29   20   18   20   22   16   16   20   21   17   11   24   16   11   18   41

If you want to do it from one season, you'd run the simulation with n_sim = 32 . . . forget about it! Here's an example:

0000 0001 0010 0011 0100 0101 0110 0111 1000 1001 1010 1011 1101 1111 
   2    2    1    2    1    1    1    4    3    1    3    4    4    3

You'll need data from multiple seasons to detect any patterns here.

3. Statistical modeling

A completely different approach is to fit a model to data. Let's assume that our helpful scrapers have done their job and supplied us with a clean dataset of all the games since 1920.

We could fit the above item-response model to the wins and losses (keeping things simple by excluding tie games).

Teams change from season to season, so I'd recommend fitting this model separately for each season, but pooling b (the home-field advantage parameter) and sigma_a (the sd of team abilities), maybe fitting a separate hierarchical model fore each decade.

But we can do better than that. Once we have the game scores, we can model them directly. Don’t model the probability of win, model the expected score differential. Something like this:

score differential ~ normal(a_i - a_j + b*home_ij, sigma_y).

The parameters a and b have slightly different interpretations now--they're on the scale of points rather than logit probabilities--but that's fine. A hierarchical model should be easy to fit. Again I'd fit a different model to each decade, or you can get more sophisticated with some sort of time series model allowing team abilities to change over the season and to have some stability between season. There's no end to the amount of modeling you can do here, if you have interest in the problem.

Once we have this model, we can simulate games and get a model-based estimate of the probability of each of the possible four-game outcomes in each season. Indeed we can do it with the actual matchups and then compare to what would happen in expectation under random matchups. The difference will give us some sense of the effect of the NFL schedule imbalances.

Once more about the z-curve method

Posted on by
36

This is Erik: I recently wrote 2 posts about my concerns about the z-curve method and some more concerns about the z-curve method. I’m sorry if this is getting repetitive, but I hope that even if you’re not interested in the z-curve per se, my criticism and especially the comments I’ve received are still interesting in light of the recent discussion about the level of rigor of meta-scientific work.

I’ll briefly summarize the z-curve method to make this post self-contained. The z-value (or z-score or z-statistic) is the estimated effect divided by the standard error, and the signal-to-noise ratio (SNR) is the true effect divided by the standard error. It is often reasonable to assume that the z-value has the normal distribution with mean SNR and variance 1.  Now suppose that we observe the z-values of a collection of studies. The z-curve method is based on the assumption that the SNRs have a discrete distribution over 0,1,2,…,6 which implies that the distribution of the z-values is a mixture of normal distribution with means 0,1,2,…,6 and unit variances. The goal is to estimate the vector of mixture weights p=(p0,p1,…,p6) and various related quantities.

The z-curve method aims to circumvent publication bias, so instead of maximizing the likelihood of the observed z-values, it maximizes the conditional likelihood of the z-values given |z|>1.96. In other words, it tries to estimate the full distribution of the z-values by using only those that exceed 1.96. The main quantity of interest is the Expected Discovery Rate (EDR) which is P(|z|>1.96) in my notation. The R function zcurve::zcurve() provides an estimate of the EDR and uses the parametric bootstrap to provide a confidence interval.

In my previous posts, I gave some examples to illustrate two concerns:

  1. The bootstrap fails as the confidence interval around p0=P(SNR=0) often collapses to the zero-length interval [0,0].
  2. z-curve can be very sensitive to slight (practically undetectable) model misspecification.

Ulrich Schimmack, who is the main author of the method, made it clear that he is unimpressed because my examples are “unrealistic” and therefore — in his opinion — irrelevant. He insists that to be realistic there should be sufficient heterogeneity among the studies, and the density of the z-values should be decreasing at 1.96. So, I did one more simulation which meets both requirements. The distribution of the z-values is the following mixture:

0.5 × N(0,1) + 0.2× N(1,1) + 0.2 × N(2,1) + 0.1× N(3,1).

The true EDR is 0.25. The coverage of the 95% (“robustified”) confidence interval across 100 simulations is 89% (CI: 81%-94%). Moreover, the estimated EDR is very biased. The average across the simulations is 0.35 (CI: 0.32-0.39). Below, I show the 100 confidence intervals. The horizontal red line is the true EDR and the blue dots are the maximum likelihood estimates (MLEs). The red confidence intervals are the cases where the confidence interval of p0 collapsed to [0,0]. It’s clear that they are an important part of the problem.

So, why is this happening? Note that the conditional likelihood of the observed z-values given |z|>1.96 has very little information about p0=P(SNR=0). That means that even if p0 is quite large, the MLE of p0 can hit the boundary, i.e. p0 is estimated at zero. If that happens, then the parametric bootstrap will create datasets without SNRs at zero. In those cases, it will be very likely that p0 will be estimated at zero. If that happens in 95% or more of the bootstrap samples, then the confidence interval collapses to [0,0].

It is actually well known in the statistical literature that the MLE can be quite biased and that the bootstrap can fail when the estimate can hit the boundary of the parameter space. See for example this paper which has a simple example.

Finally, I want to emphasize that it is not my intention to “shoot down” the z-curve method. However, I do think more work is needed before it can be used reliably. First, the problem with the collapsing confidence intervals must be fixed. Second, the limits of applicability of the method should be established and stated clearly. Third, I think z-curve is a case where the data are so weak that it’s necessary to do some regularization either with a strong prior on the mixture weights or some restrictions on them (smoothness or some shape constraint). That might be possible as several commenters on the blog seem to be very clear on what “realistic scenarios” are!

PS January 22, 2026 Frantisek Bartos uploaded zcurve version 2.4.6 to CRAN. If there are fewer than 300 z-values exceeding 1.96, the function zcurve() issues a warning:

Warning: The z-curve method is meant for large samples of test statistics. 
It might produce undercoverage and biased estimates of EDR in small sample 
sizes.

 

More concerns about the z-curve method

Posted on by
16

This is Erik: A few days ago, I wrote about my concerns about the z-curve method. I demonstrated that under certain circumstances, the coverage of the confidence interval of the expected discovery rate (EDR) is far below nominal.

Ulrich Schimmack, who is the main author of the method, posted many comments in response. I think it is fair to say that these comments are generally quite defensive and sometimes even accusatory (“Erik doesn’t want to conclude that z-curve works”). Schimmack argued that my demonstration is unrealistic and therefore not relevant. He also proposed a patch: Do not use the confidence interval for the EDR if the z-curve has an upward slope at |z|=1.96 or if the lower bound of the confidence interval of the expected replicability rate (ERR) exceeds 0.9.

Nobody likes to receive criticism, but I still find Schimmack’s reaction disappointing. If he had seriously looked at my analysis report (which I linked to in the post and shared with him and his co-authors well before the post), then he would have been able to see the main cause of the undercoverage. In 40 out of 100 simulations the confidence interval for P(SNR=0) collapses to the zero-length interval [0,0]. That’s plainly wrong. There is in fact great uncertainty about P(SNR=0).

The collapse doesn’t happen only in the “unrealistic” bimodal case which I considered in my simulation. For example, it also happens (although less frequently) when we generate samples of size n=200 from the normal distribution with mean 1 and unit variance. The interested reader can easily modify my simulation to check this and other cases.

I believe Schimmack and co-authors would do well to find out why this happens, and address the root cause of the problem. Maybe it’s just a bug in the R code. And who knows, fixing it might even make it unnecessary to “robustify” the bootstrap confidence interval of the EDR by adding ±5 percent points.

I think there’s an important general point here. If you notice a problem, or even just something weird, you should not ignore it or put a patch over it. Instead, you should try to figure out what’s causing it. In many cases, you’ll end up finding some mistake.

While I’m on the topic of the z-curve method, I would like to discuss another issue.

Confidence intervals represent sampling uncertainty – they do not take model uncertainty into account. Depending on the model, that can be a concern. The main assumption of the z-curve method is that the signal-to-noise ratio (SNR) has a discrete distribution on 0,1,2,…,6 which means that the power (probability of reaching statistical significance) has a discrete distribution on 0.05, 0.17, 0.52, 0.85, 0.98, and 1. This assumption is clearly a matter of statistical convenience. That’s fine, but it should not be expected to hold in practice. So, it’s important to see what happens when the assumption does not quite hold.

If the SNR has a discrete distribution on 1,2,…,6 then the distribution of the z-statistic is a mixture of normal distributions with means 0,1,2,…,6 and unit variances. I’ve simulated 100 samples of size n=500 of z-statistics which have a normal distribution with mean 1.5 and unit variance. Note that I’m violating the assumption of the z-curve method, but in a way that would be difficult to detect from limited data.

My report of this new simulation is here. In this case, the true EDR is 0.32. Unfortunately, the z-curve estimate of the EDR is very biased. The average of the estimates of the EDR across 100 simulations is 0.22 (CI: 0.21 to 0.23). The coverage of the “robust” 95% confidence interval is 79% (CI: 70% to 87%). Here “robust” means that ±5 percent points have been added to the bootstrap interval. The 100 (robust) confidence intervals are below; the red line indicates the true EDR.