Zombie citations: A short empirical investigation

So, we had this discussion the other day about nonexistent or irrelevant citations included in articles produced by chatbots. This is part of a larger concern about AI slop polluting the scholarly literature.

The usual way this is framed is that you look at the reference list of published article or set of articles and then you ask how many of these citations are erroneous in some way.

But maybe it’s easier to look at things the other way: Start with an article, see where it’s been cited, and follow up to see if these citations make sense.

I thought about this in the context of my classic article, published in 2010 on Arxiv and included in a book four years later, “How many zombies do you know?” Using indirect survey methods to measure alien attacks and outbreaks of the undead.

Look up this article on Google Scholar and you’ll see 6 citations. Only 6! That makes me sad. In any case, here they are:

I clicked through to read how my paper was references, and 5 of the 6 basically got it right. The sixth one messes up–it cites the zombies paper as a source for my claim that the probability is about 1 in 10 million that a vote will be decisive–I’m guessing there was some sloppy copying and pasting from my list of published papers.

So, from this short study, my estimate of the probability of a bad citation is 1/6. This was from the pre-AI era, so it’s not a “hallucination”; it’s just someone who was not careful in citation.

If the Jensen Huang / Nvidia spy scandal had happened two weeks ago . . .

If the Jensen Huang / Nvidia spy scandal had happened two weeks ago, I expect the midterm elections would’ve gone pretty much just as they actually went, but lots of pundits would’ve seized on this as the explanation for the election outcomes.

Why do I say this?

Because the result of the election was surprising. Not in its direction—everyone was expecting the Democrats to win a majority of the vote—but in the extent of the majority (in the event, they won 57% of the two-party vote for the House of Representatives) and the consequences for both houses of the legislature.

When there’s a surprise, people look for an explanation. That’s what stories are all about. As the most convenient recent politically-related news item, the Huang scandal would’ve fit the bill.

I’d still guess that, had this all happened fifteen or thirty years ago, the Huang-equivalents would’ve been Democrats. But, as it is, Huang is closely connected to an unpopular Republican president, so the story would’ve written itself. Even though I think that, had the Huang/Huang thing come up two weeks ago, the election outcome would’ve been essentially unaffected.

When is a class a “commodity” and when does it depend on the instructor?

Sometimes I get emails from students asking to work on a project with me, and I usually respond that they should take one of my classes first. If they try to work on a project without having taken a class first, I usually end up having to teach them a bunch of things. Sure, they could read my books, but I guess people don’t really like to read, also there’s something extra you get from hearing it live, doing the homework, getting feedback in class, etc.

This always puzzled me, that students would think they could jump in a project without taking any of my classes first . . . but then I realized that it’s usual to think of a math class as being pretty much the same no matter who teaches it. You take calculus or linear algebra or real analysis or whatever, and it can be taught well or poorly, there are good teachers and bad teachers, but they’re pretty much all teaching the same thing.

But statistics is different. Even if you limit yourself to a narrow topic–for example, regression analysis or survey sampling or Bayesian statistics (just to name three fairly generic-sounding classes that I’ve taught)–the course varies a huge amount from instructor to instructor. Not just teaching style (cos, yeah, my style is much different from that of every other statistics teacher I’ve seen. I’m not saying it’s better, mind you, just that it’s different. You can get a sense of it by reading Part 1 of Active Statistics), but actually what gets taught. First off, most people don’t teach regression out of Regression and Other Stories, and other regression textbooks don’t look much like ours. Second, somehow what gets covered in class isn’t just what’s in the book. This is even true with Bayesian statistics. Lots of people teach out of Bayesian Data Analysis–and I hope that, soon, lots of people will be teaching from Bayesian Workflow–but there are still so many different directions to go.

And I’m not saying I’m unique here. In general, if you want to work with statistician X, it’s probably a good idea to first take a course from X, for just about any X.

In that way, a statistics class is less like a math class or a physics class or a chemistry class and more like a history class or a literature class. Even if the books are set ahead of time, different teachers teach it so differently that they’re just different courses entirely.

But I could see how many students, whose experience is taking straight math courses in high school and college, and maybe some very empty course such as AP statistics or the college equivalent, would look at the roster of statistics classes at a university and just see a collection of topics rather than what it really is, a collection of topics interacted with instructors. And our course descriptions usually don’t do such a good job of conveying that either. Hence this post.

P.S. I’m not saying that a course being a standardized commodity is a bad thing. In a way, that what my colleagues and I have been trying to do with our textbooks, as well as with Active Statistics: we’re trying to give other people the tools to teach as we do. But, the way things stand, it still matters a lot who teaches your class–maybe not such a bad thing, as the students should get some benefit from their big tuition payments. It just seems like a secret to many students that the instructor matters, not just in how well the topic is taught but also in the content of the course.

Just Bronder it: Steve’s bronder.nyc guide to NYC off(-off)-Broadway, jazz, readings, and galleries

This post is from Bob and it is not about statistics.

Steve Bronder, our primary C++ developer for Stan these days, has built a really neat web app (nothing to download/install) consolidating all the fun stuff happening in NYC. It’s absolutely free or I wouldn’t be advertising here.

bronder.nyc: Off-Broadway and off-off-Broadway theatre, readings, jazz, and galleries in NYC right now

The web app is beautifully simple and incredibly useful. Looking for jazz shows, Mitzi and I used to surf the individual venue pages we knew, but the waters are choppy. It’s incredibly slow and frustrating given how badly most of them are laid out for information discovery. Steve’s app just lists everything out in one scrollable list that’s very easy to narrow down to a particular day or interval. The map makes it very easy to find something in walking distance. I didn’t realize just how much I needed this app until he built it and I started using it.

I would highly recommend this if you’re visiting NYC and recommend it even more if you live here. The map mode is super useful for planning things like dates. I’m not just the shill, I’m also a client. Mitzi and I just went on a double date (!) with Nawaf and Christina to the Flatiron Room in Murray Hill this past Sunday night, which was awesome. I found The Flatiron Room on the app, having never heard of the venue before. The Moses Patrou Duo cranked out some really swinging and boppy soul jazz/funk on Hammond organ and drums. Loud enough to hear, but not so loud you couldn’t have a conversation. Putting in a plug for Flatiron Room, there’s no cover, and the Japanese food was really good and reasonably priced, which is about a hundred times better than it needs to be and half the price it could be in a venue like this (this is also true of my fave sushi in NYC, Studio 151, which is above Nublu in the East Village, and has the next tier up of sushi, but still at relatively bargain prices).

P.S. I was amused by Steve’s cookies message, which is a very Steve way of stating the issue: I use Google Analytics to see which shows and theaters are popular. No personal data, no ads. [Gosh I love cookies (Accept)] [Gosh I hate cookies (Reject)].

P.P.S. I am hoping “Bronder” becomes a common word like “Zagat,” but as a verb, like “Google.” I just learned about LLM poisoning attacks, so consider this document 1 of 250.

Bayesian Workflow (my talk in NYC, Wed 18 Nov 2:30pm)

This talk is sponsored by the New York City Chapter of the American Statistical Association, and it’s at CUNY Graduate Center, 365 Fifth Avenue, room 8301 & 8304 (8th Floor).  If you want to attend (in person or remotely; I recommend in person!), just go to the form before 15 Nov and fill it out.

Also here’s a zoom link so if you’re attending remotely, maybe you don’t need to sign up.

Here’s the summary:

Bayesian statistics and statistical practice have evolved over the years, driven by advancements in theory, methods, and computational tools. This book explores the intricate workflows of applied Bayesian statistics, aiming to uncover the tacit knowledge often overlooked in published papers and textbooks. By systematizing the process of Bayesian model development, the book seeks to improve applied analyses and inspire future innovations in theory, methods, and software. It emphasizes the importance of iterative model building, model checking, computational troubleshooting, and simulated-data experimentation, offering a comprehensive perspective on statistical analysis.

Through detailed examples and practical guidance, the book bridges the gap between theory and application, empowering practitioners and researchers to navigate the complexities of Bayesian inference. It is not a checklist or cookbook but a flexible framework for understanding and resolving challenges in statistical modeling and decision making under uncertainty.

The key innovations of Bayesian Workflow include:

    • Going beyond Bayesian data analysis (model building, inference, model checking, and model expansion) to consider the larger process of statistical modeling, including comparisons of multiple models fit to a single dataset.
    • A fuller use of informative priors.  One reason we’ve moved toward stronger priors is that the replication crisis has taught us that the amount of prior information available in any given problem is often approximately the same as the information coming from an experiment (see here, for example). Informative priors also fit our increased focus on generative modeling, and we’re doing a lot more prior predictive checking to understand the implications of our models.
    • More integration between modeling, data analysis, and computing. One way to see this is that the Bayesian Workflow webpage has the code to run all our examples.
    • Simulation-based experimentation.  At its most basic, simulation-based experimentation provides a best-case baseline for statistical methods: if you can’t recover your quantities of interest with sufficient accuracy under ideal conditions (when your data are simulated from the model you’re fitting), then you know you’re in trouble. And often this is the case! Beyond that, we can simulate from one model and fit another, and see what happens. Simulation experiments aren’t always so easy to construct, as they involve specifying the entire data-generation process.

All the information on the book, including the book itself, is at the Bayesian Workflow webpage.

The non-Bayesians among you might be interested in Appendix A of the book, on Statistical and computational workflow for Bayesians and non-Bayesians, and also our short general article, Statistical Workflow, and its accompanying journal issue.

In the talk I’ll follow my usual practice and tell lots of stories.  Stories are how we learn about the world, and I’ll try to also go meta and discuss the anomalous and immutable aspects of the stories I tell.

If you can make it to the talk, please come prepared with some tough questions!

Constructing a latent-variable model to fit measurement data: A polychoric story

Following up on our recent post, Mathias Berggren writes:

I coded up a polychoric estimate for the problem you present in Constructing a latent-variable model to fit measurement data: A Stan story, part 2. I think it led to pretty decent estimates. I attach the R-code in case you want to have a look.

I hadn’t thought about it, but I guess it makes sense that this sort of model is used enough that people would have already hard-coded it!

Here’s the R code:

# Polychoric solution of latent variable stan example from Andrew

# Load required libraries
library(polycor)

# Data as a matrix
data = matrix(c(208,501,51,640,9879,1439,40,464,1432), nrow = 3, ncol = 3)

# Observed joint probabilities
Pobs = data / sum(colSums(data))

# Polychoric estimate of the latent correlation
r = polychor(data)

# Observed marginal probabilities
pX = colSums(data) / sum(colSums(data))
pY = rowSums(data) / sum(colSums(data))

# Corresponding marginal thresholds
cX = c(qnorm(pX[1]),qnorm(pX[1]+pX[2]))
cY = c(qnorm(pY[1]),qnorm(pY[1]+pY[2]))

# Estimated joint probabilities from this information
P = matrix(nrow = 3, ncol = 3)
P[1,1] = integrate(function(y) dnorm(y) * pnorm((cX[1]-r*y)/sqrt(1-r^2)), -Inf, cY[1])$value
P[1,2] = integrate(function(y) dnorm(y) * (pnorm((cX[2]-r*y)/sqrt(1-r^2)) - pnorm((cX[1]-r*y)/sqrt(1-r^2))), -Inf, cY[1])$value
P[2,1] = integrate(function(y) dnorm(y) * pnorm((cX[1]-r*y)/sqrt(1-r^2)), cY[1], cY[2])$value
P[2,2] = integrate(function(y) dnorm(y) * (pnorm((cX[2]-r*y)/sqrt(1-r^2)) - pnorm((cX[1]-r*y)/sqrt(1-r^2))), cY[1], cY[2])$value
P[1,3] = pY[1] - P[1,1] - P[1,2]
P[2,3] = pY[2] - P[2,1] - P[2,2]
P[3,1] = pX[1] - P[1,1] - P[2,1]
P[3,2] = pX[2] - P[1,2] - P[2,2]
P[3,3] = pX[3] - P[1,3] - P[2,3]

# Compare estimated and observed probabilities
P
Pobs
Pobs - P

# Find least squares estimates to transform thresholds to 1.5 and 3.5
b1X = sd(c(1.5, 3.5)) / sd(cX)
b0X = mean(c(1.5, 3.5)) - b1X * mean(cX)

b1Y = sd(c(1.5, 3.5)) / sd(cY)
b0Y = mean(c(1.5, 3.5)) - b1Y * mean(cY)

# Check  so we get the correct values
b0X + b1X * cX
b0Y + b1Y * cY

# Least squares estimates
# These show how we should transform the standard normal marginal 
# distributions above to a scale with thresholds 1.5 and 3.5.
# Thus, they also indicate the updated mean and standard deviation of 
# each variable (as they transform from one with mean = 0 and sd = 1).
# If one wants one common mean, then one might take some weighted mean 
# of the intercepts
c(b0X,b1X)
c(b0Y,b1Y)

# Discussion
# If there are more thresholds then the least squares solutions may not be 
# exact, and we may find that the thresholds cannot be at 1.5, 2.5, etc.
# Same way, one might find that the intercept cannot be the same.

The challenge of truly probabilistic programming

Many years ago, Jouni and I wrote an unpublished paper called Fully Bayesian Computing:

A fully Bayesian computing environment calls for the possibility of defining vector and array objects that may contain both random and deterministic quantities, and syntax rules that allow treating these objects much like any variables or numeric arrays. Working within the statistical package R, we introduce a new object-oriented framework based on a new random variable data type that is implicitly represented by simulations.

A few years later we published an updated version of the paper under the title, Manipulating and summarizing posterior simulations using random variable objects. Here’s the background:

In the Bayesian paradigm, unknown quantities have probability distributions and are thus random variables. Observed values are just realizations of random variables, and constants may be thought of as random variables with point mass distributions. In mathematical notation, we deal with objects that are random variables, but in practice these objects are approximated by vectors of numbers, that is, simulations. Consequently, when programming for manipulating simulations of unknown quantities, we must write code to manipulate arrays of numeric constants.

And here’s the challenge:

Arrays of simulations are cumbersome objects to work with. Functions that work with vectors will not in general work with matrices, so special versions of the functions need to be written to accommodate matrices of simulations as arguments. For example, a scalar-valued random variable becomes a vector of simulations, and a random vector becomes a matrix of simulations. . . . The usual approach now is to use loops and vector-matrix computations to manipulate the matrix of posterior simulations. This approach is general but awkward and far from transparent . . .

As part of this effort, Jouni wrote a package in R called “rv” that allows direct manipulation of random variable objects, so for example if alpha and beta are random variables, each a vector of length J, and gamma is a scalar, then you can compute something like delta = (alpha + gamma)/beta, which would be a vector of length J. Internally, alpha, beta, and delta would be stored as S x K arrays, and gamma would be stored as an array of length S, where S is the number of simulation draws in your simulation frame, but the user would be able to work with these without having to keep track of that indexing; it would all be implicit.

The rv package was just a demonstration, but then more recently Paul Bürkner, Jonah Gabry, Matthew Kay, and Aki Vehtari implemented the idea in the “posterior” package in R. They call it rvar: The Random Variable Datatype. I still haven’t integrated rvar() into my Bayesian workflow, but I’d like to do it.

The point I want to make here is that this random variable datatype moves us closer to what I’d call “truly probabilistic programming.”

Here’s what I mean by that.

Stan (and Bugs, PyMC, etc.) are called “probabilistic programming languages,” characterized by Wikipedia as “based on the declarative specification of probabilistic models, for which inference is performed automatically.” Fair enough. In Stan, the variables defined in the parameters block, transformed parameters block, model block, and generated parameters block represent random variables with a joint density that is specified by the Stan program.

But, although this is called a probabilistic program as the term is used, it doesn’t fully feel like probabilistic programming in the sense of allowing me to grab random variables and work with them directly. I can do all that within the Stan program in the generated quantities block, but at a certain remove. Stan and other probabilistic programming languages are set up to run in batch mode: you specify a model and then the inference engine runs and spits out posterior simulations (or some other distributional summary).

To me, a truly probabilistic programming environment would allow me to interactively work with random variables, and to work with them as random variables, not as arrays of simulation draws. I want to be able to multiply two matrices without having to loop through the simulations. Part of this is about clean code, but more generally it’s about workflow. Any mathematical or graphical operations we can do on numbers, we should be able to do with random variables, with the simulations all hiding under the hood. Kind of by analogy to how we can do floating-point math on the computer without (usually) having to worry about how all these digits are being stored.

This will require tools such as rvar() and also some experience with the workflow. The payoff is . . . you could call it “better postprocessing of Bayesian simulations,” but I think that would understate the potential benefits of being able to work directl with uncertainty without that extra array dimension getting in your face at every step.

“You’re kind of a hero at my company”

This post is by Phil Price, not Andrew. 

The work I’m proudest of is not necessarily the work that is most useful or best known or that others consider the best. A few years ago I learned that a model I threw together in about 15 minutes, and whose simplicity I had been quietly embarrassed by for years, was used by a successful company.  But I just learned that the model is much more widely used and appreciated than I thought. And I’ve just lost my two steady consulting clients, who were together responsible for about 60% of my income each year, so now is the time to do some self-promotion. And I think there are broader lessons here about the virtues of simplicity, and how we sometimes aren’t good judges of our own work, and maybe some other themes besides. So here’s the story.


About three years ago, a friend and colleague named Steve Greenberg died of cancer, only 64 years old. Steve was a very hands-on kind of guy (as you can learn from his obituary, ) with a deep practical understanding of HVAC systems, which made him very valuable in our quasi-academic setting where a lot of us understood the theories but wouldn’t know the first thing about how to actually install or manage the equipment. Steve was also a committed bike commuter. At his “celebration of life” I wore a bike jacket and a hard hat to recognize those aspects of Steve. Of course my hard hat has my name on it, not Steve’s name.

In the drinks and snacks portion of the party, after the speeches and slideshows, a stranger came up to me and said “Phil Price! You’re kind of a hero at my company.” This was, to say the least, surprising. OK, I’ve done some consulting work for various companies, but none of them would consider me a “hero,” not even remotely close, so what could this possibly be about? But it turns out they weren’t talking about my consulting work at all, but about some work I had done at LBNL years ago: inventing the “Time-of-week and Temperature” (TOWT) model to predict building electric load. “That model is the basis for most of what my company does to measure and verify energy-saving measures in buildings.” That’s not an exact quote, I’m sure, but it’s how I remember it.

Here is an explanation of the model that is clearer and more compact than anything I ever wrote about it. Or, if you want to see more detail and some extensions of the model and how to use it for various purposes, here’s a video of talk I gave back in 2015 (11 years ago!) as I was leaving LBNL to become a consultant. (For me, the video doesn’t display right in Safari but works fine in Firefox, not sure what’s up with that).

The TOWT model is useful in spite of its simplicity, or perhaps even because of its simplicity. It’s very easy to explain, it’s pretty robust, it’s very easy to code up and to fit, it is easy to extend in various ways (e.g. by including additional explanatory variables). But I never bring it up, it has barely been mentioned on my CV (although I’m fixing that), and I haven’t previously bragged about / emphasized my role as creator of the TOWT model, because I’m kinda sheepish about it. As you can see from the web page above or the video, the model is extremely simple: there’s a time-of-week indicator variable, and a piecewise-continuous linear function of temperature…and that’s it, at least in the simplest version of the model. I literally came up with it in an afternoon; actually less than an afternoon, it was more like 15 minutes from conception to implementation. It is very simple and I’ve worried that if I use it as an example of what I’ve done, people will think that level of sophistication is all I can do!

But maybe that’s a silly thing to worry about. People had been analyzing the same kinds of data for decades and in all that time nobody had come up with this model (as far as I know!), simple though it is. And jeez, I had no idea it was so widely used, it’s not just that one company. This is what a search turns up, for tools and methods that rely on the TOWT model:

  •  RM&V (Reference Measurement & Verification Tool): An open-source M&V reference tool released by LBNL that implements piecewise linear regression using TOWT and machine learning methods for multi-site facility analysis.
  • NMECR: An R-based package by kW Engineering implementing peer-reviewed Normalized Metered Energy Consumption (NMEC) algorithms, featuring an expanded TOWT model that supports operational schedules.
  • OpenEEmeter / OpenDSM: An open-source analytics toolkit (formerly OpenEEmeter, now operating via OpenDSM) that incorporates the CalTRACK hourly methodology, which is directly adapted from LBNL’s TOWT model for demand response and energy efficiency baselines.
  • Energy Twin (ET Analytics / ET Core): A commercial energy analytics platform and Python/R library package (ecpEnergyTwinAnalytics) that uses custom TOWT-based implementations to derive operational KPIs and anomaly detection across commercial portfolios.
  • Power TakeOff (Verify™): An automated M&V 2.0 platform that uses TOWT-based modeling to isolate weather versus schedule loads for utility energy efficiency programs.
  • CalTRACK Hourly Methodology: A standardized protocol for calculating normalized metered energy savings using an adapted hourly TOWT architecture for commercial and residential buildings.
  • NMEC (Normalized Metered Energy Consumption): A protocol used by utilities (such as in California) where the TOWT model serves as a benchmark standard to quantify weather-normalized energy savings.
  • TOWST (Time-Of-Week, Solar and Temperature): An extension of the standard TOWT model that integrates physical constraints like solar gain and building orientation via iterative clustering to improve prediction accuracy.

So, OK, I’m in. That’s my model. I’ll try not to be sheepish about it anymore.

As far as the general lessons to be learned, I think there are a few:

  1. Simplicity can be good. So what if the TOWT model was originally conceived and implemented in a few minutes. Indeed, that’s a good thing! If a simple model does the job, use the simple model.
  2. Don’t hide your light under a bushel basket. If you’ve done something that people find useful, it doesn’t matter if you don’t think it’s your best work or if it doesn’t demonstrate the full range of your skills. This shamelessly self-promotional post was prompted by comments from people who looked at my LinkedIn profile and the letters I’ve sent to potential clients and told me I wasn’t nearly clear enough about what I’ve accomplished. In trying not to come across as bragging, I was selling myself short.
  3. You often don’t know which of your work matters. I had no idea TOWT was in widespread use until a stranger told me at a memorial gathering. I would have guessed that some of my more elaborate work was my most influential. It wasn’t… as far as I know (but how would I know?). It’s worth writing up even the simple stuff (and publishing it, or anyway putting out there somehow).
  4. “It only took 15 minutes” leaves out the work before it. The model took 15 minutes to conceive and implement because I’d spent a couple of months looking at building energy data and knew what mattered and what didn’t, and of course I was proficient in linear regression and other stats methods. Maybe I under-valued my work because I forgot about this fact.

This post is by Phil

Survey Statistics: Sampling from a Bayesian Menu

A week ago Andrew wrote about Alan Zaslavsky, who recently passed. Alan’s obituary has information about an upcoming celebration of life.

Alan Mark Zaslavsky Profile Photo

Andrew highlighted Alan’s work on the census. In Sampling from a Bayesian Menu, Alan wrote about objections not only to model-based adjustments for undercount, but the Supreme Court’s ruling in 1999 against even purely “design-based” estimates using sampling for nonresponse followup.

Alan was on the panel who wrote “Sampling for Nonresponse Follow-up: Achieving Adequate Precision at Acceptable Cost”. They explained that after initial response to the census, the bureau conducts nonresponse follow-up on remaining housing units. In the 1990 census, this was massively delayed, costing both money and accuracy. They recommended that the 2000 census should conduct nonresponse follow-up on a random sample of nonresponding housing units for unbiased estimation of the total number of people in all nonresponding housing units. This is what Alan described in Sampling from a Bayesian Menu as the “design-based” approach that was rejected by the Supreme Court in 1999.

Despite this bad (old) news, the paper is fun to read, and will make you hungry.

Constructing a latent-variable model to fit measurement data: A Stan story, part 2

Following up on our post from the other day, we were trying to come up with parameters for the following measurement-error model:

– People have latent a distribution of attractiveness:
latent_attractiveness ~ normal(0, sigma_latent)
– Attractiveness for each person is measured twice, with different mean levels and different levels of added variance:
At time t, observed_attractiveness ~ normal(mu_obs_t + latent_attractiveness, sigma_obs_t), for t = 1, 2
– Attractiveness at each time is then recorded as a discrete variable, rounded to the nearest integer and constrained to fall between 1 and 5: That is, if observed attractiveness is less than 1.5, it is recorded as 1, if observed attractiveness is between 1.5 and 2.5, it is recorded as 2, etc.
– The only information easily accessible is recorded attractiveness binned into three categories: 1 or 2, 3 or 4, and 5. This summary is available for each person at both time points.

Here are the data, from the first two waves of the Add Health survey:

So, here’s the deal. We want to choose five parameters (sigma_latent, mu_obs_1, mu_obs_2, sigma_obs_1, sigma_obs_2) to fit 9 data points. Actually 8 data points because what we’re fitting are the proportions that sum to 1.

This shouldn’t be so hard, right?

Actually, it is kind of hard! At least, it’s hard to get a fit that gets real close to the data table.

1. Trial and error to find parameter values that match the data

I started with trial and error, first computing the mean and sd of the attractiveness ratings for each wave of the data (approximating the three binned levels as 1.5, 3.5, and 5 on the 1-to-5 scale), then setting mu_obs_1 and mu_obs_2 to the mean values and partitioning the variance so that sd_latent and sd_obs are roughly equal. In R:

discretize <- function(a) {
  ifelse(a < 2.5, 1, ifelse(a < 4.5, 2, 3))
}
N_table <- rbind(c(208, 640, 40), c(501, 9879, 464), c(51, 1439, 1432))
N_data <- sum(N_table)
long_data <- array(NA, c(N_data, 2))
long_data[,1] <- rep(c(1,2,3, 1,2,3, 1,2,3), as.vector(N_table))
long_data[,2] <- rep(c(1,1,1, 2,2,2, 3,3,3), as.vector(N_table))
long_data_1to5 <- ifelse(long_data==1, 1.5, ifelse(long_data==2, 3.5, ifelse(long_data==3, 5, NA)))
means <- apply(long_data_1to5, 2, mean)
sds <- apply(long_data_1to5, 2, sd)
mu_obs <- means
sigma_latent <- sqrt(0.5*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)

The result:

> print(round(mu_obs, 2))
[1] 3.68 3.59
> print(round(sigma_latent, 2))
[1] 0.54
> print(round(sigma_obs, 2))
[1] 0.61 0.45

I then estimate the probabilities in the 9 cells by simulating a million pairs of observations from this model, first simulating latent attractiveness, then the two observed attractiveness measures, then rounding:

N <- 1e6
y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)

And here's the result:

> print(table(y_discrete[,1], y_discrete[,2]))
         1      2      3
  1  18715  54792    393
  2  41167 674479  55053
  3    653 110776  43972

OK, it's hard to compare this to the data table because the numbers are so big, so let's rescale so that it sums to the number of people in the data:

> print(table(y_discrete[,1], y_discrete[,2])*N_data/N)
              1           2           3
  1  274.249610  802.921968    5.759022
  2  603.261218 9883.815266  806.746662
  3    9.569062 1623.311504  644.365688

This is hard to read so I'll round it:

> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
       1    2    3
  1  274  803    6
  2  603 9884  807
  3   10 1623  644

And, for convenience, let's look again at the data:

> print(N_table)
     [,1] [,2] [,3]
[1,]  208  640   40
[2,]  501 9879  464
[3,]   51 1439 1432

Comparing the cells [1,3] and [3,1] in the simulated table to the observed data, we see that our simulations don't have enough variation between the two measurements. So let's rejigger by decreasing sigma_latent and increasing sigma_obs:

sigma_latent <- sqrt(0.25*mean(sds^2))
sigma_obs <- sqrt(sds^2 - sigma_latent^2)

Simulate again with these new parameter values:

y_latent <- rnorm(N, 0, sigma_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- rnorm(N, mu_obs[j] + y_latent, sigma_obs[j])
}
y_discrete <- discretize(y_obs)

And here's what we get:

> print(round(table(y_discrete[,1], y_discrete[,2])*N_data/N))
   
       1    2    3
  1  145  896   37
  2  687 9592 1011
  3   54 1827  404

Now the [1,3] and [3,1] entries are basically OK, but [1,1] and [3,3] are too small.

I tried some more tinkering with the five parameters but couldn't get the table to line up. And then I thought: Hey, why not cut out the middleman and just solve the optimization problem directly?

2. Stan program: first try

It goes like this:

data: the observed table, N_table
parameters: the five parameters, sigma_latent, mu_obs, sigma_obs
model: compute table(y_discrete) and compare to N_table

In this model, table(y_discrete)*N_data/N represents E(N_table|parameters), and we were computing this using simulation . . . but you can't simulate random data in the Stan model block! Or we could compute the expectations directly, but these are two-dimensional normal integrals, and there's no formula for that.

So I did a hack. You can't do simulations in the Stan model block, but you can simulate random numbers offline and pass them into the Stan model. Here's our Stan program, table_match.stan:

functions {
  int discretize (real y) {
    if (y < 2.5)
      return 1;
    else if (y > 2.5 && y < 4.5)
      return 2;
    else
      return 3;
  }
}
data {
  array[3, 3] int N_table;
  int N;
}
transformed data {
  int N_data = 0;
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_data += N_table[j1,j2];
    }
  }
  vector[N] z_latent;
  matrix[N, 2] z_obs;
  for (n in 1:N){
    z_latent[n] = std_normal_rng();
    for (j in 1:2) {
      z_obs[n, j] = std_normal_rng();
    }
  }
  real a = 1.*N_data/N;
}
parameters {
  vector[2] mu_obs;
  real s_latent;
  vector[2] s_obs;
}
transformed parameters {
  matrix[3,3] N_sim = rep_matrix(0, 3, 3);
  for (n in 1:N){
    real y_latent = s_latent * z_latent[n];
    array[2] int y_discrete;
    for (j in 1:2){
      real y_obs = y_latent + mu_obs[j] + s_obs[j] * z_obs[n,j];
      y_discrete[j] = discretize(y_obs);
    }
    N_sim[y_discrete[1], y_discrete[2]] += 1;
  }
  N_sim = a*N_sim;
}
model {
  mu_obs ~ normal(3,3);
  s_latent ~ exponential(1);
  s_obs ~ exponential(1);
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_table[j1, j2] ~ poisson(N_sim[j1, j2]);
    }
  }
}

At the very end we're using the Poisson likelihood to compare the observed data to the model. And we include weak priors on the model parameters to keep everything stable.

We run the program from R:

library("cmdstanr")
table_match <- cmdstan_model("table_match.stan")
table_data <- list(N_table=N_table, N=1000)
fit <- table_match$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit, max_rows=30)

Aaaand . . . it doesn't work:

   variable     mean   median       sd      mad       q5      q95  rhat ess_bulk ess_tail
 lp__       77677.57 77880.76 10570.22 13247.66 63077.85 91926.34  3.66        4        5
 mu_obs[1]      0.24     0.10     0.61     0.51    -0.42     1.22  3.94        4       11
 mu_obs[2]      0.47     0.40     1.64     2.24    -1.47     2.54  3.70        4       11
 s_latent       2.78     1.82     2.52     1.37     0.49     7.00  3.42        4       12
 s_obs[1]       1.84     1.60     1.18     1.08     0.46     3.68  3.47        4       13
 s_obs[2]       1.65     1.70     1.36     1.93     0.18     3.03  3.43        4       11
 N_sim[1,1]  8290.04  8865.67  1305.70   640.92  6096.06  9319.94  5.54        4       NA
 N_sim[2,1]   890.16   593.49   669.06   358.48   351.70  2007.60 32.33        4       NA
 N_sim[3,1]   468.94    80.60   711.48    65.18    14.65  1699.86 13.53        4        4
 N_sim[1,2]  1759.81  1472.73  1368.41  1292.70   219.81  3883.31  5.95        4        4
 N_sim[2,2]   852.42   622.80   655.06   369.34   219.81  1948.98  9.66        4       NA
 N_sim[3,2]   322.47   307.73   251.19   336.75    29.31   644.78  8.22        4        4
 N_sim[1,3]   626.65   256.44   779.88   260.71    14.65  1963.64 18.31        4       NA
 N_sim[2,3]   408.84   300.41   339.12   271.58    58.62   952.51 16.39        4       NA
 N_sim[3,3]  1034.67   381.00  1331.88   304.16    29.31  3326.46 16.98        4       NA

What about Stan in optimize setting?

fit_optimize <- table_match$optimize(data=table_data)

Nope:

Initial log joint probability = 69172.8 
    Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
       1       69172.8             0       9.10807    0.001749       0.001       45    
Chain 1 Optimization terminated with error: 
Chain 1   Line search failed to achieve a sufficient decrease, no more progress can be made
Warning: Fitting finished unexpectedly! Use the $output() method for more information.

Damn.

3. Stan program: revision

What went wrong? The above code is mathematically correct, but the discreteness of all these simulations leads to difficulty with the Hamiltonian Monte Carlo. The problem is that, due to the use of fixed error terms and a finite N, continuous changes in the parameters result in jumps in the expected data and thus jumps in the likelihood.

But then I realized something: I can't compute two-dimensional normal integrals in closed form, but one-dimensional normal integrals are trivial. So, instead of passing in the vectors of random numbers z_latent, z_obs[,1], and z_obs[,2], I can just pass in z_latent, and then we can use the normal cdf function in Stan to integrate over z_obs[,1] and z_obs[,2]. This offers two advantages: first, by integrating over two of the three dimensions analytically rather than by simulation, we're computing the cell probabilities more accurately and we don't need such a large N; second, the target function will now be a smoothly continuous function of the parameters, so HMC should work better.

Here's our revised Stan program, table_match_2.stan:

data {
  array[3, 3] int N_table;
  int N;
}
transformed data {
  int N_data = 0;
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      N_data += N_table[j1,j2];
    }
  }
  vector[N] z_latent;
  for (n in 1:N){
    z_latent[n] = std_normal_rng();
  }
  real a = 1.*N_data/N;
}
parameters {
  vector[2] mu_obs;
  real s_latent;
  vector[2] s_obs;
  real sigma;
}
transformed parameters {
  matrix[3,3] N_sim = rep_matrix(0, 3, 3);
  for (n in 1:N){
    real y_latent = s_latent * z_latent[n];
    matrix[3,2] p;
    for (j in 1:2) {
      p[1,j] = normal_cdf(2.5 | y_latent + mu_obs[j], s_obs[j]);
      p[2,j] = normal_cdf(4.5 | y_latent + mu_obs[j], s_obs[j]) - p[1,j];
      p[3,j] = 1 - (p[1,j] + p[2,j]);
    }
    N_sim += p[,1] * p[,2]';
  }
  N_sim = a*N_sim;
}
model {
  mu_obs ~ normal(3,3);
  s_latent ~ exponential(1);
  s_obs ~ exponential(1);
  sigma ~ exponential(1);
  for (j1 in 1:3) {
    for (j2 in 1:3) {
      sqrt(N_table[j1, j2]) ~ normal(sqrt(N_sim[j1, j2]), sigma);
    }
  }
}

You can see from the last line of code that I also changed the likelihood. The Poisson likelihood is appropriate if the model is correct, but it won't be correct; by fitting on the square root scale but including a free parameter, sigma, we're allowing for overdispersion.

We call it from R:

fit_2 <- table_match_2$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_2, max_rows=30)

Aaahhhh, much better:

   variable    mean  median     sd    mad      q5      q95 rhat ess_bulk ess_tail
 lp__        -25.68  -25.29   2.22   2.08  -29.93   -22.80 1.00     1009     1687
 mu_obs[1]     3.79    3.79   0.08   0.07    3.67     3.92 1.00     2090     2047
 mu_obs[2]     3.66    3.66   0.07   0.06    3.56     3.78 1.00     1762     2092
 s_latent      0.68    0.68   0.06   0.05    0.59     0.77 1.00     2562     1803
 s_obs[1]      0.51    0.51   0.08   0.07    0.39     0.64 1.00     2562     1894
 s_obs[2]      0.27    0.28   0.11   0.10    0.08     0.44 1.00     1791     1139
 sigma         4.22    4.07   1.01   0.93    2.84     6.14 1.00     1511     2172
 N_sim[1,1]  338.69  326.65 114.39 103.68  173.55   539.20 1.00     2433     2009
 N_sim[2,1]  434.37  418.15 164.48 152.14  199.93   722.66 1.00     2775     2271
 N_sim[3,1]    1.63    0.52   3.59   0.70    0.01     6.71 1.00     2207     2454
 N_sim[1,2]  588.17  570.14 193.09 180.59  304.30   926.40 1.00     2580     2098
 N_sim[2,2] 9734.95 9745.35 520.62 488.97 8858.11 10563.91 1.00     3395     2533
 N_sim[3,2] 1637.06 1621.71 322.01 300.94 1134.19  2184.29 1.00     3159     2436
 N_sim[1,3]    0.66    0.13   1.92   0.18    0.00     2.98 1.00     2309     2220
 N_sim[2,3]  595.11  575.31 193.51 175.52  319.79   954.80 1.00     2496     2634
 N_sim[3,3] 1323.36 1298.78 300.80 278.37  863.93  1847.16 1.00     1967     1755

Although, if you look at the expected values in the table (N_sim), it's not so great! The [1,3] and [3,1] values are 0.66 and 1.63, much lower than the data values of 40 and 51.

I think the problem is with the Poisson or square-root likelihood, that it allows relatively large error for small-count cells.

4. Stan program: another revision

So I'll try a new model, this time with a lognormal likelihood. This will go into the file table_match_3.stan. I'll just show you the one line in the model block because that's the only part of the program that's changed:

      log(N_table[j1, j2]) ~ normal(log(N_sim[j1, j2]), sigma);

We call it from R:

fit_3 <- table_match_3$sample(data=table_data, parallel_chains=4, refresh=0)
print(fit_3, max_rows=30)

Which yields:

   variable    mean  median      sd     mad      q5     q95 rhat ess_bulk ess_tail
 lp__         -6.98   -6.42    3.35    3.34  -13.32   -2.51 1.01      339      559
 mu_obs[1]     4.12    4.00    0.59    0.37    3.48    5.25 1.00      819      721
 mu_obs[2]     3.86    3.74    0.52    0.22    3.38    4.89 1.01      761      600
 s_latent      0.82    0.69    0.48    0.22    0.40    1.69 1.01      434      306
 s_obs[1]      0.98    0.95    0.37    0.28    0.46    1.58 1.00     1040      744
 s_obs[2]      0.55    0.51    0.34    0.23    0.09    1.16 1.01      777      833
 sigma         0.84    0.74    0.38    0.30    0.42    1.56 1.01      448      810
 N_sim[1,1]  488.28  286.72  716.36  190.52   91.35 1536.07 1.00     1071      888
 N_sim[2,1]  666.97  552.96  480.97  294.04  220.03 1474.32 1.00     1772     1300
 N_sim[3,1]  115.33   64.10  308.25   42.04   19.43  302.69 1.01     1794     1046
 N_sim[1,2] 1043.18  944.89  607.16  461.48  313.54 2133.92 1.00     1109     1153
 N_sim[2,2] 5820.30 6119.15 2055.18 2028.03 1832.99 8729.56 1.01      432      610
 N_sim[3,2] 2887.52 2773.01 1308.13 1259.89  941.96 5226.20 1.00     1142     1137
 N_sim[1,3]   69.16   36.91  145.99   24.02   12.50  203.99 1.00     1730     1098
 N_sim[2,3] 1139.58  969.93  709.41  479.43  390.12 2416.68 1.00     1478     1611
 N_sim[3,3] 2423.68 1726.31 2045.37 1221.81  495.16 7118.74 1.01      526      635

Or, in optimize setting:

fit_optimize <- table_match$optimize(data=table_data)
print(fit_3_optimize, max_rows=30)
   variable estimate
 lp__           1.62
 mu_obs[1]      3.86
 mu_obs[2]      3.68
 s_latent       0.54
 s_obs[1]       0.82
 s_obs[2]       0.52
 sigma          0.38
 N_sim[1,1]   209.05
 N_sim[2,1]   549.51
 N_sim[3,1]    49.04
 N_sim[1,2]   931.32
 N_sim[2,2]  8175.10
 N_sim[3,2]  2736.90
 N_sim[1,3]    30.44
 N_sim[2,3]   984.07
 N_sim[3,3]   988.56

Not bad! Not perfect---you can see that the [3,3] cell is much too big. But it'll do.

5. Conclusion

The model doesn't quite fit! Maybe it could make sense to try longer-tailed distributions such as the t with low degrees of freedom. But for our original purpose, which is to understand the impact of measurement error on a nonlinear relationship, this will be just fine. Also we've learned something about statistical computing.

Can peer review survive unrestricted AI use by reviewers?

This is Jessica. Reviews from ICLR (one of the three biggest ML conferences) came out last Friday, and a large number of them appear to have been primarily AI-generated, at least if we believe what authors are reporting on X. I believe it. On my own submissions, about half had Pangram scores over 50% (with some papers getting majority reviews with 90+% Pangram scores). There are other signs as well, like appendix mention reportedly being nearly twice as common in reviews compared to last year. 

Was this rampant AI use in review surprising? The official review policy for ICLR 2027 permitted “limited and responsible use of AI tools to assist human writing of peer reviews and meta-reviews,” but said it was “mandatory to disclose use of AI tools,” with the review form soliciting this information and asking that if you use AI you also report your original self-written paper assessment and any LM interactions. 

It’s a little mind-boggling to me that anyone would still think asking people to disclose AI use is a viable policy. By now it should be common knowledge that when you say something like “mandatory to disclose use of AI tools” that will be interpreted as “suggested to disclose,” and no one will do it.  

It’s also a bit mind boggling that no one anticipated this mess. At ICML 2026, many reviewers who deliberately selected into a pool that would not be allowed to use AI were caught using AI by watermarking. Turns out that when you pump high numbers of papers at each reviewer their good intentions start to crumble. And ICLR received a record number of abstract submissions – over 60,000, suggesting total valid submissions probably surpassed last year’s total near 20,000. 

When reviewers feel comfortable telling authors “I have used AI tool to do the review since I do not have enough time”* something is broken. Good for that reviewer for being honest, but who is sending the impression that an AI review you produced but can’t understand is a viable way to fulfill your reviewer agreement? 

This situation makes a few things clear. One is that it’s irresponsible to give reviewers high review loads but not restrict how reviewers can use AI, at least in a field like ML where many reviewers are inexperienced and required to do multiple reciprocal reviews. 

The second thing is that it’s irresponsible not to prioritize preserving review over hurting authors’ feelings. ICLR did introduce a cap of 20 submissions per author, plus a limit of one paper per author if none of the authors qualified for reviewing. But it’s too conservative. Compare to the policy adopted by STOC, a CS theory conference, which this year capped submissions at five per author, dropped anonymity, required authors to submit a video in which they explain their paper, and reserved the right to desk reject for poor exposition. We can debate whether these are the right steps, but at least someone on that organizing committee had some balls. 

Why letting external reviewers use AI freely is a bad idea 

At this point, every author is (or should be) using AI to check their own work if they care about quality. Frontier models can find errors that you and your co-workers miss. It can suggest limitations that reviewers might bring up. I don’t send out papers anymore if I haven’t run them through Refine.ink.

So what is the point of reviewers, with less expertise on what the authors did than they have (and in ML, we’re often talking junior PhD students), meat-proxying to send them more AI-generated feedback? Seemingly none.

Consider what we lose if we ban AI in peer review. The main risk would seem to be to peer review’s supposed function as a soundness filter. LLMs enhance a reviewer’s ability to comb through 50+ pages of paper and appendix to verify that what is presented is accurate. But even if the authors are lazy and miss some problems themselves, letting individual reviewers feed papers to LLMs using basic prompts to “write a review for [conference X]” is far from the best solution. Strong LLMs too can miss major errors when not directed what to look for (see e.g., examples from Nihar Shah’s work).

Error identification and basic soundness checks are much better addressed by the venue itself, using validated AI tools. Centralize checks like whether key results can be reproduced from the provided data/code, verify whether equations, proofs, or calculated statistics are correct, detect leakage or violated assumptions, make sure key claims can be traced to evidence, and so on. There’s precedent for all of this in existing research, see Table 1 here or arguments others have made for using AI in review for such purposes here and here.

Some of these automated checks can be used to determine if a paper should be desk rejected. Some will be more controversial (and they should be, as I discuss in more detail here, because there are no hard-and-fast filters to separate good from bad science). Anything contestable can be surfaced for reviewers. 

Naturally this can shift what cues get attention during review. When we give reviewers a shared summary representation of the paper, we correlate errors relative to what we’d get if they relied on their private signals alone. But these are problems that machine learning researchers can define and study! See e.g., our work formalizing what it means to have a good LLM-implemented checklist (i.e., use review corpora to learn a set of checks that reduce the combination of bias and variance under LLM-assisted review relative to the variance you’d have from human-only review). We can also educate reviewers on what their role is–-figure out what matters that can’t be verified directly across the paper (significance, novelty, etc.) and decide how to weight the surfaced checks. 

But can’t good reviewers write better reviews with AI? 

One counterargument to banning reviewers from using AI at their discretion is that we are crippling review quality. It’s easier to get more information from a paper if you can query an LLM about it while you read. A reviewer may reach better insights if they can rely on LLMs to help them run simulations or produce counterexamples to understand the method or theory. It’s not that all AI use is bad. It depends on whether the critiques or suggestions that AI surfaces are within the convex hull of the reviewers’ expertise. If they are, the value of their review might be amplified by AI use. But when they’re not, we get reviewers confidently reporting claims from an LLM that they cannot themselves verify, which the authors may already have decided were not good advice in their own interactions with AI. This is where review becomes pointless. 

To retain the value of AI for engaging with papers, venues may also want to provide access to LLMs for interacting with the paper. This makes it possible to preserve some provenance in where claims in reviews are coming from, which authors and ACs can use to appraise reviews as needed. It becomes possible for venues to run checks that reviewers are contributing some kind of unique information to the process, rather than parroting directly. If reviewers are adding minimal information, then those reviews should be removed. 

Ultimately, we need reviewers to be using their own judgment to weigh in on the significance and importance and novelty of ideas. Otherwise fields may drift in directions that are misaligned with what the humans working in those fields value. The problem with unchecked AI use is that even when people might want to stay in control of their ultimate judgments, they are not always good at recognizing when they are taking on biases reflected by LLM output. There was a study out of Anthropic recently where they endowed LLMs with preferences that made them lean toward particular recommendations, then had people use them to do a task. The majority of the time, people adopted the preference that was endowed to their LLM, and many of them also attributed the recommendations to themselves. 

Why are we still protecting authors’ rights to submit slop?

Reviewers turning to AI en masse is a symptom of conferences being hesitant about desk reject policies. For example, at the top ML conferences, one can still submit a mostly AI-produced paper and expect a team of human reviewers to provide feedback. I’ve heard many students who are reviewing for ICLR this year complain about how many papers in their review piles were hard-to-read AI slop.

Organizers add some minimal friction for authors, and expect the reviewers to suck it up and deal with the remaining onslaught. But when you don’t prioritize reviewers’ ability to actually do their job by requiring more authorial accountability, you hasten the “review death spiral”, where as review gets noisier, authors are further incentivized to submit low quality papers, in the hopes that they get lucky. 

There are better ways to deal with this problem than the ultra-conservative don’t-hurt-any-feelings approaches we’ve been seeing. We should not be seeing teams of undergrads submitting 9 papers at a time only to withdraw them once they see the initial reviews, as reported by one AC. For example, start desk rejecting papers that are flagged by Pangram as mostly AI-written, or that are very hard to read. If an author can’t take the time to signal that there is a human on the other end who is accountable for the work, do we really owe them 3+ other humans’ time providing feedback? 

I really hope that policies tighten. As imperfect as peer review has always been, removing people from it entirely is not going to lead us to better science.

*From a screenshot of a review shared on X before the author took it down.

Elizabeth Holmes, Jensen Huang, etc.: If all this had happened fifteen to thirty years ago, I’m guessing that all these people would’ve been Democrats. What changed?

We recently saw You Can See Everything—yes, I highly recommend it, as with The Odyssey, it’s as good as everybody says, it lives up to the hype—and then there was the recent news about tech executive Jensen Huang. Put these two stories together and you get the above charming bit of incoherence. “Bring back the gallows,” indeed. Billy Evans seems to be discriminating in his thoughts on the recommended severity of punishment for white-collar crime.

In the movie, Billy Evans, the husband of disgraced medtech fraudster Elizabeth Holmes, says that he’s a fan of the late conservative firebrand Charlie Kirk, and, as you can see below, the couple are all-in on Donald Trump. Huang too is closely linked to the president, who in turn is a promoter of Wild West ventures such as crypto and betting markets.

While watching You Can See Everything and reading the news, I had this thought: If all this had happened fifteen to thirty years ago, I’m guessing that all these people would’ve been Democrats. Big tech money is right-coded today but was left-coded in the (Bill) Clinton and Obama eras.

The question is, why has this happened? The trend is particularly notable when juxtaposed with the general shift of college-educated voters away from the Republican party during the past two decades.

One difference, perhaps, is that information technology is more closely tied to the military than it was twenty years ago. Back then, the defense-contractors and tech were connected, but I think less intertwined. So it could make sense for Apple, Google, Microsoft, etc., to side against George W. Bush and the Iraq War, while current AI companies don’t want to go against a militarily-active government. Taking over Greenland is more exciting to today’s tech investors; twenty years ago they were more interested in opening up consumer markets around the world, and an unpopular foreign invasions were not a useful step in that direction.

Another take on it is that business investors are inclined to go with the winner. In 2008, this was Obama; in 2024, it was Trump. And I do think this is part of the story but not the whole thing. In 2008, Obama’s team had the high-tech luster and the Republicans were the grumpy outsiders. Since 2016, it’s been the opposite. Even setting aside Jensen Huang and Elizabeth Holmes, not to mention Elon Musk and that space mirrors guy, there’s a general sense that conservatives are now the tech cheerleaders and liberals are the new Debbie Downers.

There’s a story to explain any given individual case, or, more generally, you can say that U.S.-style conservatives support something closer to an unregulated free market, while U.S.-style liberals are more concerned about the balance of power between business and labor, so in that sense this all makes sense. One theory that’s been floating around is that tech executives have moved to the right in reaction to efforts by the left-wingers in their technical workforce to reduce executives’ freedom of action. Go back twenty or thirty years, and the Democratic party was led by figures such as Bill Clinton, Al Gore, and Barack Obama, all of whom were friendly not just to big business but to the tech industry in particular. This brings us closer to political scientist Thomas Ferguson’s view of national politics as reflecting divisions within business executives in competing industries. Also, tech lords were rich a generation ago—remember all those jokes about Bill Gates being rich—but not nearly as rich in aggregate as they are now, so that has to be part of the story too.

It’s also interesting to compare how different people react to right-coded business fraud vs. left-coded academic fraud.

“Causal Effects Versus Causal Mechanisms: Two Traditions With Different Requirements and Contributions Towards Causal Understanding”

Following up on our post from last year on the importance of mechanistic models, Jim Grace points to his above-titled paper in Ecology Letters. From the abstract:

When methods are described as ‘causal,’ emphasis is increasingly placed on statistical techniques for isolating associations so as to quantify causal effects. In contrast, natural scientists have historically approached the pursuit of causal knowledge through the investigation of mechanisms that interconnect the components of systems. . . . the case is made that an explicit assessment of existing mechanistic knowledge should be an initial step in causal investigations.

Grace also has this webpage.

I’m posting this because I continue to be interested in the connections between mechanistic and statistical models of causality, and I continue to think that the fields of statistics, economics, and policy analysis took a wrong turn by associating “causal inference” with black-box statistical analysis. I understand the appeal of the idea–indeed, the very first applied statistics class I ever took, 40 years ago, was Rubin’s course on causal inference–but in many ways it’s been a very long detour that only in recent years have we been able to take in a more useful direction.

2026 compared to other dramatic midterm election reversals of the past half-century

As I said a few days ago, you don’t need me anymore to do your fancy post-election data analyses. So instead I can offer speculation (Could the 2026 election results lead to a decrease in political polarization?) and meta-commentary (NYT columnist offers epically stupid take on 2026 midterm elections).

But we can also learn from some really simple analyses.

Comparing 2026 to previous midterm swings

How big of a deal was the 2026 national election?

To look at consequences, you can count the number of seats won by the two parties. To look at the aggregate opinions of voters, it makes sense to look at the national popular vote for the House of Representatives, where the Democrats received 56.2% of the two-party vote, a swing of 7.5 percentage points from their House vote share in 2024.

Is 56.2% with a swing of 7.5pp a little or a lot? One way to judge this is to compare to the midterm stunners of recent decades:

2010 is the election that put the brakes on the Obama government: the Tea Party movement convinced the voters that government expansion was going too far. In that election, the Republicans won 53.5% of the two-party vote, a 9.0 percentage point swing from 2008.

2006 was the year that George W. Bush’s post-9/11 high finally wore off. The Democrats’ share of the two-party vote in that election: 54.1%, a 5.5pp swing from the previous election.

1994 was the year that Newt Gingrich and the Republicans offered a clear alternative to the liberal policies of the Clinton administration, taking control of Congress with 53.5% of the two-party vote, a 6.1pp swing in their favor. 53.5% is lower than the other numbers I’ll be showing today, but it arguably understates the Republicans’ accomplishment in 1994, in that this was their first stint in the majority in forty years.

1982 isn’t remembered much now, because the Republicans roared back two years later, but, with the worst economic crisis since the Great Depression, and President Reagan having very low job approval, the Democrats did very well at the ballot box, winning 56.0% of the two-party House vote, a 4.1pp swing from the previous election.

1974 was the post-Watergate election, bringing many fresh faces into Congress. The Democrats capitalized on economic stagnation and discontent with the criminal doings of the Nixon administration to gather 56.9% of the two-party vote, a 4.2pp swing.

In summary, the Democrats’ performance in 2026—a 7.5 percentage point swing yielding 56.2% of the two-party vote—was indeed exceptional, at or above the level of previous repudiations of the incumbent party during the midterm elections.

This can be understood in various ways.

One way that 2026 was different is that this time the incumbent party was really the incumbent party. Unlike in those previous campaigns, this time the Republicans had unified control of the federal government—the legislative, executive, and judicial branches—and the incumbent party was itself highly unified (as can be seen, for example, in the Senate confirmation of some extreme cabinet officers). So, voters who wanted party balancing had a clear choice in a way that wasn’t so obvious in previous years where there was divided government.

The other way is to compare to historical baselines. The Democrats’ victories in 1974 and 1982 cemented their status as the default majority party, but that’s not as impressive as the Republicans’ feat in 1994, and then in 2010, of moving from the minority to gain firm control of Congress. The Democrats achieved this in 2026.

What comes next?

In most of the above cases, the midterm winner did not follow up with a decisive win at the next presidential election. 2006, followed by Obama’s victory two years later, was the only exception. Consider the other cases:
After 1974, the Republican party was left for dead: it had been the majority party for awhile and then was wiped out in the congressional elections. But, two years later, the not-particularly-popular Gerald Ford came within 2 percentage points of being elected president.
After the repudiation of Republicans in 1982, the economy came back to life—as the president quipped, now that things are going better, “they don’t call it Reaganomics anymore”—and the presidential election was a landslide.
And the Republicans’ congressional victories in 1994 and 2010 were followed by Clinton and Obama doing just fine in their subsequent reelection campaigns.

That said, the Republicans’ successes in 1994 and 2010 represent one possible template: the triumphant return of a party that had been given up for dead. In both cases, the victorious congressional party did not follow up with a presidential election win, but the party remained strong enough to determine much of the political agenda going forward.

The take-home message based on this history: the Democrats should value their congressional win for itself and not assume they’re on a glide path for 2028. Conversely, the Republicans should figure out how best to work with the Democratic majority without thinking of themselves as being about to lose all power.

Nationalization of politics

Along with partisan polarization has come a nationalization of the vote, which has led to two things. First, the vote swings in presidential elections have been much smaller. After their losses in 1982, the Republicans rebounded a bit in 1984, upping their House vote share to 47.4%, a gain of 3.4 percentage points. In that same year, Reagan was reelected with a two-party vote share of 59.2%!

That sort of thing doesn’t happen anymore, in two ways. First, with partisan polarization, you wouldn’t expect either party’s presidential candidate to come close to 60% of the vote. Second, congressional elections are nationalized, and you wouldn’t expect to see a stunning 12 percentage point gap between a party’s congressional and presidential vote share.

Constructing a latent-variable model to fit measurement data: A Stan story, part 1

This simulation is a little different than what came before. For my earlier post I simulated hypothetical future data from the beauty-and-sex-ratio study. For my new post I’m simulating data from a hypothetical strong relationship between two variables. Here I’m entirely focused on the question of what happens to a nonlinear (but monotonic) pattern when independent measurement error is added to the predictor. To make the patterns clear, I will assume a large effect; indeed, I’m only concerned here with the effect size, not its variation.

The one thing I’ll preserve from the previous simulation is the distribution of attractiveness measurements in the sample: 2%, 5%, 45%, 37%, and 11% in the five attractiveness categories.

I want to set up a model where there’s an assumed true nonlinear relationship between some outcome of interest and a continuous underlying latent attractiveness variable.

The nonlinear model isn’t hard to come up with. The hard thing is creating a distribution for the latent attractiveness, under the constraint that it be roughly consistent with the data.

I’ll assume a symmetric model, a normal distribution on the measurement scale. This might seem like a strong assumption but it seems reasonable enough given the data, which have an average that’s between 3 and 4. The tails are asymmetric but that’s just because the scale has 5 as an upper bound.

My first step is to compute the mean and sd of the observed data:

n <- round(2792*c(.02,.05,.45,.37,.11))
mean_ratings <- sum(n*(1:5))/sum(n)
sd_ratings <- sqrt(sum(n*((1:5) - mean_ratings)^2)/sum(n))

These come to 3.5 and 0.83. Assuming random zero-mean measurement error, the latent attractiveness distribution should have the same mean but a lower variance. We should subtract 1/12 to account for the rounding (as this is approximately equivalent to adding an independent uniform(-0.5,0.5) random variable, which has a variance of 1/12) and then more to account for measurement error.

But how much to subtract? For this, we need to make some assumption about variation in the attractiveness measurement. Attractiveness seems pretty subjective to me (also it can vary a lot over time, as we know from those high-school yearbook photos of people who later became movie stars). I gotta assume something here, so let's assume that the standard deviation of the measurements is 0.5. I'm picking this because 0.5 is a large value, but it's not as large as the sd of 0.83 in the data. Indeed, 0.5^2 + 1/12 = 0.33 and 0.83^2 = 0.69, so I'm assuming that about half of the variance in the data is coming from measurement error.

But is 0.5 a reasonable value? To answer that question, it helps to have data.

The data at hand come from wave 3 of the Add Health study. I already downloaded those data awhile ago in preparing my earlier analysis, but I've had difficulty working with the numbers. I can't remember why, but there's no single dataset and no single codebook. It's kind of a mess. I spent a few minutes going through to see if there were any repeated measurements of attractiveness there, but I wasn't sure where to look.

So I googled *reliability of add health attractiveness measure*. I was hoping to find the names of the relevant variables in the survey, and that didn't come up in the search. What did appear was a bunch of research papers on the topic. I looked at a few of them and then came across this one: "Ugly criminals," by Naci Mocan and Erdal Tekin. I'm not so interested in the topic of the paper---I'll leave it to others to decide if its conclusions are reasonably supported by its data---; what's relevant here is that they discuss the Add Health attractiveness measure. Which is indeed how it came up in my search.

They looked at data from respondents whose attractiveness had been recorded in two waves. That's exactly what we need for a reliability analysis:

There's only one mild annoyance here: they pooled categories 1 and 2, and they pooled categories 3 and 4. So we don't have the full 5 x 5 table, only a 3 x 3 collapsing of it.

But we can deal with that. We'll just simulate and see what happens.

First we write a couple of simple functions:

round_and_bound <- function (a, lower=1, upper=5){
  pmin(upper, pmax(lower, round(a)))
}

clumped <- function(a) {
  ifelse(a==1 | a==2, 1, ifelse(a==3 | a==4, 2, ifelse(a==5, 3, NA)))
}

We set up the parameters for the latent data:

mu_obs <- c(3.87, 3.69)
s_latent <- 0.53
s_obs <- c(0.82, 0.52) 

Then we simulate two independent measurements on each of 14654 people---that's the number of cases in the above table:

y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}

And then we look at the two margins and the table showing the joint distribution of the two measurements:

print(table(clumped(observed_1)))
print(table(clumped(observed_2)))
print(table(clumped(observed_1), clumped(observed_2)))

I also looked at the unclumped data but I won't bother showing you that here.

Here are the results:

> print(table(clumped(observed_1)))
    1     2     3 
  668 11186  2800 

> print(table(clumped(observed_2)))
    1     2     3 
  723 11221  2710 

> print(table(clumped(observed_1), clumped(observed_2)))
       1    2    3
  1  227  437    4
  2  490 9292 1404
  3    6 1492 1302

These don't quite match the table above. The simulations have not enough "unattractive" people and too many "very attractive." Also they two waves in the simulations are essentially the same (by construction, as we simulated them from the same distribution), but the real data show a few more unattractive people and a lot more very attractive people in wave 1.

There would be a way to set up my simulation to alter its parameters and get the closest possible match to the observed two-way table---I guess it wouldn't be hard to program this in Stan, actually! But here I took the lazy way out and just tinkered with the numbers to get things to work. It took about five tries of playing with the inputs, and eventually I converged on this:

sd_measurement <- 0.5
mu_latent <- 3.8
sd_latent <- 0.67
latent_beauty <- rnorm(n_sim, mu_latent, sd_latent)
observed_1 <- round_and_bound(rnorm(n_sim, latent_beauty, sd_measurement))
observed_2 <- round_and_bound(rnorm(n_sim, latent_beauty - 0.1, sd_measurement*0.6))

As you can see, I hacked the simulations to give the wave 2 measurement a lower mean and lower variance.

Here's the result:

> print(table(clumped(observed_1)))
    1     2     3 
  845 10839  2970 

> print(table(clumped(observed_2)))
    1     2     3 
  731 11852  2071

> print(table(clumped(observed_1), clumped(observed_2)))
       1    2    3
  1  330  514    1
  2  400 9727  712
  3    1 1611 1358

OK, now I have the margins pretty good, but the joint table doesn't have enough people whose measurements shifted between the lowest and highest categories.

I played around some more with the parameters of the model but I couldn't get the counts to line up with the data.

So I bit the bullet and did what I should've done from the start, which was to write a Stan program to find a best fit. That was its own struggle---as is often the case, I learn a lot from these little examples---so I'll devote a separate post to it.

Here's what I ended up using:

mu_obs <- c(3.87, 3.69)
s_latent <- 0.53
s_obs <- c(0.82, 0.52) 
y_latent <- rnorm(N, 0, s_latent)
y_obs <- array(NA, c(N, 2))
for (j in 1:2) {
  y_obs[,j] <- round_and_bound(rnorm(N, mu_obs[j] + y_latent, s_obs[j]))
}

And here's the result, which is still not perfect (as will be discussed in my followup post) but is enough for my goals here:

> for (j in 1:2) print(table(clumped(y_obs[,j])))
   1     2    3 
 1225 9684 3745  

   1     2    3 
 815 11857 1982 

> print(table(clumped(y_obs[,1]), clumped(y_obs[,2])))
       1    2    3
  1  208  990   27
  2  563 8141  980
  3   44 2726  975

Which is far from a perfect match to the actual table of responses from the two waves:

> print(N_table)
     [,1] [,2] [,3]
[1,]  208  640   40
[2,]  501 9879  464
[3,]   51 1439 1432

But it's close enough, as we're just trying to make a general point about attenuation of effects under measurement error.

P.S. For the next story, see Constructing a latent-variable model to fit measurement data: A Stan story, part 2.

Could the 2026 election results lead to a decrease in political polarization?

I have some hope, for the following reasons:

1. The incoming Democrats in Congress include some centrists and some on the left. Mix them in with the ones who are already there, and we can expect divisions on general ideological grounds and on specific issues such as Israel/Palestine, data centers, etc.

2. There have always been congressional Republicans who detest Trump, but they’ve shown close to zero defections from the president’s policies. But now that they’re the minority, the forces keeping them in line will be weaker. Also they are getting the point that general elections as well as primaries are threats to them.

My hope is that there will be some important issues that split within each party, leading to cross-party coalitions. I’m not saying this would be a cure for polarization—back in 1993, the vote on NAFTA split the two parties, and that did nothing to stop the Newt Gingrich-led polarization that took over soon after—but maybe it would be a step forward.

To put it another way, I think of there being two components to partisan polarization. One component is a natural consequence of connections between issues: there’s a logic to the left and right representing different coalitions and being represented by two political parties. The other component is not so good, it’s the lockstep conformity that leads to candidates taking extreme positions to appeal to primary election voters, presidents crafting policies to satisfy moderates just enough to get the threshold of required votes, and party-line votes to confirm unqualified cabinet officials. Maybe some open disputes on high-profile issues will make pure partisan voting less of an expected thing, so that legislators can vote against their party leadership without it being such a big deal, which in turn could motivate moderation among their leaders. Which does not mean that non-centrist policies on the left or right could not be passed, just that they’d require coalition building of some sort..

Political polarization has been steadily increasing for the past thirty-five years, so it seems like it would be a mistake to bet against more of the same. But trends can occasionally reverse; nothing lasts forever. Back when my colleagues and I were starting to do research in political science, the incumbency advantage had been steadily increasing for forty years . . . but then the trend reversed.

Also, partisan polarization is a funny thing: some people are strongly motivated to vote for Republicans because they hate the Democrats and what they stand for; others are strongly motivated to vote for Democrats because they hate the Republicans etc. There aren’t too many voters who actually like either party, at least not on a national level. To borrow a physical analogy, you could say that political polarization is maintained by tension, not compression.

P.S. Somewhat relatedly, NYT columnist offers epically stupid take on 2026 midterm elections.

P.P.S. I think I blew it by not using a clickbait title such as “Polarization is so over.”

“So in that sense the silly statistics in the paper is no big deal, it serves as a sort of rhetorical flourish?”

Richard Feldman writes:

I’m a long-time reader and sometimes commenter on the blog. I’ve learned a lot – and been entertained – by many of the posts.

A paper was just published in Science on a theme that is potentially of interest to your blog readership, and that has the potential to enter the news, given its political nature. It’s not my research area, but the statistical practice to which I object crops up in my discipline – ecology – quite often. In fact, I’m reviewing a paper right now that does something similar.

The paper is “Partisan disparities in the funding of science in the United States.”

What has me question my sanity as someone trying to engage in good statistical practice is that the authors ran a regression with the covariates of interest (party control of the presidency, house, and senate) and many “control” variables and then a simple model with the covariates of interest and just one “control” variable. (A whole other issue is using the term control to describe the additional covariates and to refer to party control…in the same figure!). In the main figure in the text, they present both results – from the simple and full model – but clearly emphasize the results from the simple model. They justify their emphasis on the simpler model with what sounds ludicrous to me: “Collinear controls may absorb relevant variance, and the descriptive meaning of a regression coefficient becomes increasingly opaque as more variables are included in the model.”

I call the above the “have your cake and eat it too” approach to interpreting regressions. You propose that there are a series of factors to consider when measuring the effect of a variable of interest, but then interpret the effect as if those other factors weren’t important. I think there tends to be a default sentiment that coefficients in multiple regression are difficult to interpret – in the paper I’m reviewing, this is what the authors say – because of a lack of understanding that it is the effect when the other variables are at a pre-specified value (usually zero or their mean, but now with R packages like marginaleffects it is easy to specify those values). Am I right to have an unease about the approach of running multiple and simple regressions at the same time, but interpreting just the simple regression?

In addition to the above, the paper analyses an untold number of models – “alternative specifications” – and uses an NHST approach without considering the inflation of Type I error.

I agree that this is a pretty empty paper, published in a prestigious journal just because it’s on a topic of interest to the editors and also it’s wrapped in a clean-looking package. In addition to the flaws noted above, there’s also the problem that they don’t start their analysis earlier, a period during which the Democrats were in continuous control of Congress for several decades. On the other hand, their general point seems correct, that until this year science funding has had steady bipartisan support. So in that sense the silly statistics in the paper is no big deal, it serves as a sort of rhetorical flourish?

Analysis of tonight’s election results (and why I don’t need to do it anymore)

Back in the day, I would do a post-election post, summarizing information from the election outcomes and exit polls. For example, after the 2008 election:

As with previous Republican candidates, McCain did better among the rich than the poor:

outcome1.png

The gap between young and old has increased—a lot:

ages.png

But there was no massive turnout among young voters. According to the exit polls, 18% of the voters this time were under 30, as compared to 17% of voters in 2004. (By comparison, 22% of voting-age Americans are under 30.)

The red/blue map was not redrawn; it was more of a national partisan swing. See this state-by-state scatterplot of Obama vote in 2008 vs. Kerry vote in 2004:

2004_2008_actual.png

The standard deviation of the state swings (excluding D.C. and the unusual case of Hawaii) was 3.3%. That is, after accounting for the national swing in Obama’s favor, most of the states were within 3% of where they were, compared to their relative positions in 2004.

By comparison, here’s the 2000/2004 graph:

2000_2004.png

The standard deviation of these state swings was 2.4%. This was even less variation–2004 was basically a replay of 2000—still, the relative state swings of 3.3% in 2008 were not large by historical standards.

Again, Obama didn’t redraw the map; he shifted the map over in his favor. (Or, to put it more precisely, the economy shifted the map over in the Democrats’ favor and Obama took advantage of this.)

Here’s the map showing where Obama and McCain did better or worse than expected based on 2004:

2004_2008_map_actual.png

Finally, how did the pre-election polls do? Here’s the scatterplot:

2008_2008.png

I see some systematic patterns: Obama underperformed where the polls had him way down, and he outperformed where the polls had him up.

But I don’t need to do this sort of thing anymore.

Why not? Because the news media do it well all by themselves now. The New York Times page is good, also freelance sites by analysts such as former Fivethirtyeight.com honchos Nate Silver and Elliott Morris. There’s been an impressive improvement in political analysis in the past 25 years.

Back in the day, the election-night and post-election summaries in the news media were frustratingly crude. As the votes came in, the TV networks would just report raw totals based on X% of the vote in state Y, without even trying to do any analytics, even anything as simple as tallying the swing from those precincts compared to the previous election. And then the next day you just get a list of the winners along with some exit poll results.

What’s changed?

First, exit polling has become more difficult—all those mail-in votes—so there’s more of a focus on the actual election results. Second, outsiders such as Nate have upped the game, leading other news organizations to follow. Third, data processing is better so that now sites can immediately make this sort of display:

Things didn’t use to be this way. Back in 2000, I was in a meeting with the news media exit poll consortium and I was trying to persuade them to look at swings and to do MRP to help with their noisy data . . . and they didn’t listen to me. They were a happy monopoly and didn’t need no fancy analytics. But look at the news media now. I have the impression that campaigns are much more sophisticated in their analyses too.

I’d like to say that the media and political campaigns have caught up with academic political science. But it would be more accurate to say that they’ve surpassed us. I still think academic political science has useful things to contribute: new methods such as MRPW, innovative data collection and presentation from CES and others, and also we have the space to try to systematically understand trends such as trickle-down polarization. As well as all sorts of work out there that I haven’t even heard of. But these analyses I used to do, looking at swings, going back and forth between the local story and the big picture? The news media can do it better. I guess that’s how Bill James must feel when he sees modern sports analytics: he has his insights, but the stuff he did that was difficult and innovative back in the 1980s is now out there, and at a much more sophisticated level.

P.S. I’m hoping in a few months to get a detailed report from Yair and his team, reporting not just what happened in the vote but what was going on underneath. Kind of like their reports from 2018, 2020, 2022 and 2024. Each of these came out in the May following the election, so stay tuned.

Survey Statistics: Times/Siena Polls and Times Voter Trends

Yesterday Andrew wrote about Alan Zaslavsky, who recently passed.

This is typical of his work in that it was innovative, with applied importance, and understated.

As Fan Li quoted, Alan was “the best applied statistician on earth”. I’ll write more Alan stories soon. He was my PhD advisor and hero.

Since it is election day, let’s talk about methods to analyze polls, where Alan’s legacy can’t be overstated. Fan Li mentioned work with Alan on balancing covariates, cited by Huling & Mak (2024) in their work on energy balancing weights, which are used by the Times/Siena Poll. This Survey Statistics series has discussed a few Times/Siena Poll decisions over the last year.

This week, in “The Final Polls Show a Blue Surge. Can We Trust Them?” Nate Cohn compares how the 2026 polls would look without these big changes, and sees they would be on average more than 3 points higher for Democrats.

Last week Andrew wrote about MRP at the NYT. So cool ! Their methodology article says:

The estimates do not report a margin of sampling error because they are not based on a sample of the population. The estimates are derived from multiple models, which are each subject to different kinds of error.

Yes, there are many kinds of error to consider, in “Margin of Error” we discussed how to get an MRP margin of error assuming no nonresponse bias after the MRP adjustments, no measurement error, and no model misspecification. In “Total Margin of Error” we saw that in election polls from 1998 to 2014, Shirani-Mehr et al. found that total survey error is about twice the reported margins of error. Even if you have margins of error to report, visualizing them is challenging (see Gelman and Price 1999 “All maps of parameter estimates are misleading”).

I was curious about the relationship between their MRP methods and their weighting methods. (See “MRPW”, “struggles with equivalent weights”, continued struggles, and “equivalent models, equivalent weights (locally)”, and “structured MRP to smooth survey weights”.) The only mention I found of their weights is:

The polls were joined to the given year’s voter file and reweighted to match the voter file by L2’s modeled race, party, age and interactions thereof, with the original survey weight as the base weight. The adjusted survey weight was used in all subsequent race models.

So they calibrate the base weights to match the voter file, see Lumley 2010 Chapter 7.

It’s been a few months since I’ve blogged about the logit shift, so it is time. (See: “2 flavors of calibration”, 2nd helpings (multinomial outcomes), 3rd helpings (multivariate outcomes), 4th helpings (more on multivariate outcomes), “logit shift and raking”.) The NYT uses a version of this shift:

The resulting individual-level probabilities were adjusted in multiple steps to match the actual election result by precinct. First, the probability that a voter would support a minor party candidate was adjusted to match the minor party proportion of the vote. Second, the probability that a voter would support a major party candidate was adjusted to match a candidate’s major party vote share in the precinct. Finally, the probability that a voter would support the Democratic or Republican candidate was calculated by multiplying the adjusted Democratic and Republican share of the major party vote by the adjusted major party vote share of the precinct.

Presumably they’re doing something like:

  1. E(shift_a(P(major party | X))) = P(major party) solve for a.
  2. Then either:
    1. E(shift_b(P(Dem | major party, X))) = P(Dem | major party) solve for b. Or
    2. E(shift_a(P(major party | X)) * shift_b(P(Dem | major party, X))) = P(Dem) solve for b.

2A isn’t quite right: it implies E(shift_a(P(major party | X))) * E(shift_b(P(Dem | major party, X))) = P(Dem). But we want 2B to hold, and the two are equivalent only if the individual-level probabilities for major party and for Dem given major party are uncorrelated.

It’s hard to write prose to describe something technical. The NYT level of methods transparency is great !