Survey Statistics: how to measure hospital quality ?

Last week we chatted about hopes for a new edition of Lumley 2010.

Preview

Alan Zaslavsky recommends this textbook in his course. After taking the course in 2011, I asked Alan to advise my PhD. Alan’s work extends across most topics in Lumley’s textbook and beyond. One focus is statistical issues in measuring the quality of a hospital: how can we compare hospitals who see different distributions of patients (i.e. casemix) ? Zaslavsky 2001 explains a few approaches, including direct standardization.

Let:

  • Y = complications
  • X = sickness of the patient
  • R = hospital, can be A or B
  • Y(A) = complications if treated by hospital A
  • Y(B) = complications if treated by hospital B

To compare hospitals A and B, we want to compare Y(A) versus Y(B). We cannot do this for one patient, because at least one of Y(A) or Y(B) is not observed. But we can compare averages. We can compare E(E(Y(A) | X, R = A)) versus E(E(Y(B) | X, R = B)) where the outer expectations are with respect to some standard casemix of X. Direct standardization estimates the inner expectations with the mean in each X stratum. As Zaslavsky 2001 points out, this needs modification if there are no patients in a stratum, or when X is continuous. Instead, the inner expectations can be estimated with regression models.

See Figure 2 from Zaslavsky 2001:

Lumley 2010 points out that direct standardization is mathematically equivalent to poststratification:

 

Survey Statistics: textbooks adding non-probability sampling chapters

Last week I asked for your favorite survey statistics textbooks. Thank you for many fun recommendations here and on Bluesky ! I am compiling a reading list for after the election. Here’s some not-blue-sky on the AT in TN:

Ben Schneider and Winston Lin share my love for Lumley 2010. I agree that it’s fun and a second edition would rock.

Preview

Some recommended books added sections on non-probability sampling in their new editions:

Lumley 2010 already includes material applicable to non-probability sampling, see e.g. “2 flavors of calibration”. If there is a second edition maybe it will follow the books above and add a chapter about non-probability sampling ?

(See “probability samples vs epsem samples vs SRS samples” for more about the term “non-probability sampling”.)

Survey Statistics: 2 different uses of simulated data

Last month Andrew wrote:

I use simulated data to understand how a model works; it can be viewed as a form of mathematical analysis, a way of working out the implications of a set of assumptions.

I don’t use simulated data to learn about the external world. For that, I need to take measurements of reality.

Let’s start with the latter, which concerns using chatbot-simulated survey responses.

Simulated data to learn about the external world:

Andrew wrote: “the chatbot is trained on the past, but the usual reason for surveys is to learn about changes”. Pew Research agrees: “Survey-taking is best left to humans”. And a year ago we saw that Thomas Lumley agreed: “It [Claude or ChatGPT] will probably do this better than I could, but it’s not magic.”

That said, Andrew suggested using the simulated response as an additional adjustment variable. Jessica Hullman wrote last year about a suggestion from Broska, Howes, and van Loon (2025) that uses LLM predictions as auxiliary data in something like a survey regression estimator. Andrew of course suggests Multilevel Regression and Poststratification (MRP) to incorporate auxiliary data like LLM predictions.

Now back to the first use of simulated data.

Simulated data to understand how a model works:

Andrew has written lots about this. I teach it in my course with Dr Arjun Potter at the Nelson Mandela African Institution of Science and Technology (NM-AIST). Our example is Potter et al. (2026), who study the effects of drought, fire, and herbivory on growth of various acacia tree species. In our course, we simulated data ahead of real data collection to understand the model, e.g. how block effects change the precision of the fire effect. Students understood that this exercise does not replace real data collection. We are hoping to teach these ideas again in the new year (funding dependent). For the past three iterations, see “Hand-drawn Statistical Workflow at Nelson Mandela” and “Survey Statistics: connections to experimental design”.

Survey Statistics: ANOVA follow-up

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).

Survey Statistics: Fat Bear Week 2026

Happy Fat Bear Week to all who celebrate. In 2024 I made a cartoon called Basu’s Bears, adapted from Basu’s (1971) elephants example, a lesson on the use of auxiliary information in survey statistics. For Fat Bear Week 2025, I wrote accompanying R code.

Gaurav’s post points out:

What makes the estimator wobble is not HT’s algebra but the sampling plan.
…
If a proxy x can be measured …. use probability‑proportional‑to‑size (PPS)
…
A complementary route…is to model y | x

As we saw in “equivalent models, equivalent weights (locally)”, Little 2004 linked these complementary routes, showing under which outcome model Horvitz-Thompson does well:

When is outcome model (8) correct ? When bear size y is proportional to probability we choose to measure them pi plus noise. If we follow Gaurav’s good advice and use PPS sampling, pi is proportional to auxiliary data x. So the outcome model is correct when bear size y is proportional to auxiliary data x plus noise. If x is bears before feasting, this means their expansion ratios thru salmon consumption are constant plus noise.

But what goes wrong in Basu’s (1971) example ? Instead of PPS, the compromise design is to set pi = 99/100 for Sambo and = 1/4900 for all other bears. So the outcome model is only correct when Sambo eats all other bears except for crumbs. If we can anticipate this, then HT even in this sampling plan will do just fine.

 

Survey Statistics: ANOVA

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.

Survey Statistics: Pew Research finds “No Easy Fix for Bogus Respondents in Online Opt-In Polls”

I just read the Pew Research Center’s new report “No Easy Fix for Bogus Respondents in Online Opt-In Polls”. It begins:

One of the most urgent problems in online opt-in polling is bogus (or fraudulent) respondents. These are survey-takers who make no effort to answer questions truthfully and instead are just looking to finish surveys quickly and collect rewards.

To address this, they try 3 screening methods: 1) trap questions, 2) prescreening, and 3) voterfile matching. There is no direct way to know how many bogus respondents are removed with each method. So instead they compare them with 3 data quality measures: 1) “yea-saying”, 2) quality of open-end responses, and 3) response order effects.

The “yea-saying” data quality measure is the % who say “yes” to at least 10 of 15 questions. Without screening, this was 7%. At first this seemed fine to me and I was confused why this is even a data quality measure. But a probability sample estimates only 1% of people say “yes” to this many of these questions. Screening with trap questions helped the opt-in sample match the probability sample. One of their trap questions was asking if folks use a made-up social media platform called Fizzypress.

They made a separate set of calibration weights (see “3 flavors of survey weights”) for the unscreened sample and for the 3 screening method subsamples, calibrating to ACS along these dimensions:

How does screening for bogus respondents affect results for Harris-vs-Trump 2024 vote choice ?

Say we want E(Y), Harris vote choice in the population. But we only observe the opt-in sample, so we have nonresponse error (beyond Pew’s weighting adjustment): E(Y | R_opt_in = 1) – E(Y). Also, we only observe Y* != Y due to measurement error from bogus responses.

Say everyone we screen-in has no measurement error. Then the error with screening is nonresponse: E(Y | R_opt_in = 1, screened_in) – E(Y). This may be bigger than nonresponse error without screening: E(Y | R_opt_in = 1) – E(Y), but this isn’t knowable. We only know overall error including measurement error: E(Y* | R_opt_in = 1) – E(Y). Pew found that overall error was less without screening, screened-in folks overrepresented Harris.

How do opt-in samples compare to probability samples ?

Pew writes that there is “no way for potential bad actors to self-select into probability-based samples.” Of course, a randomly selected person can still give bogus responses (responding with Y* instead of their true Y). But the rate of measurement error in probability samples would match the population: P[Y != Y* | R_probability_sample = 1] = P[Y != Y*]. Whereas the rate of measurement error in opt-in samples is presumably higher due to folks selecting into the survey to earn rewards: P[Y != Y* | R_opt_in = 1] > P[Y != Y*]. In other words, the rate of measurement error is subject to nonresponse error, a mixing of the two sides of the Groves et al. figure below.

Image

In this Survey Statistics series we’ve seen measurement error in an adjustment variable X, e.g. recalled vote: see “is a mismeasured X better than none at all ?”, “more adventures in mismeasured X”, “more on recalled vote”, “it is (still) the people”. We saw that the NYT used to drop adjustment for recalled vote due to its measurement error, which could worsen nonresponse error. Pew is studying the effect of screening out people with a lot of measurement error, which could worsen nonresponse error.

In this series we’ve also seen measurement error in the outcome variable Y, e.g. disease status, see “wanting workflow”. This example modeled the measurement error, rather than screen out mismeasurements.

Also relevant is Andrew’s post “She wants to know what are best practices on flagging bad responses and cleaning survey data and detecting bad responses. Any suggestions from the tidyverse or crunch.io?”.

Survey Statistics: logit shift and raking

We’ve been discussing how to use population margin information about poststratification variables (see Modeling Complex Contingency Tables and SynthMargins & Bayes-Raking). I want to connect two terms for this we’ve mentioned a lot in this series: the logit shift and raking.

In June 2025 we discussed 2 flavors of calibration, including the logit shift: calibrate p(z | x, survey) to match known aggregates p(z).

In June 2025 we discussed 3 flavors of survey weights, including calibrated weights: find weights w such that weighted survey E(wx | survey) matches known totals E(x). When we have at least 2 variables, say x and z, we can calibrate their cross-classification, or one margin at a time. The latter is called raking.

We can think of raking as finding weights w(x,z) or equivalently finding population probabilities p^(x,z) = w(x,z) p(x,z | survey) that match known totals. In the logit shift we want to find the population conditional probability p^(z | x) that matches a known total. These seem really similar ! Let’s compare.

Here’s my setup:

  • 2 variables: binary z, categorical x
  • Have p(x,z | survey), and p(x), p(z) population margins
  • Want p^(x,z) = argmin_q KL(q, seed) that matches known margins
    For raking: seed = p(x,z | survey)
    For logit shift: seed = p(z | x, survey) p(x)

So here’s some homework: show that Lagrange multipliers gives p^(x,z) = seed(x,z) exp(- 1 – d_x – d_z). And so the solution p^ will have the same odds ratio as the seed (which is the same for raking as the logit shift).

We can substitute these solutions into the constraints, and for either seed we get:

sum_x p(x) logit^-1( logit p(z=1 | x, survey) – d_z ) = p(z=1)

This is one equation with one unknown shift d_z. In other words, logit shift and raking give the same solution. This can be solved with a root-find, as Andrew commented that Lei et al 2017 did in this code.

More generally, if z is multinomial, we can’t get such a closed form equation with one unknown, so we can use Iterative Proportional Fitting (IPF), as we discussed in 2nd helpings of the logit shift.

Survey Statistics: more on SynthMargins and Bayes-Raking

Last week we discussed SynthMargins, a method from the poster Modeling Complex Contingency Tables that uses partial information (margins) about poststratification variables. On theme for this series (“it is the people”), the authors commented thanking Andrew for introducing them, folks from different fields with a shared goal: Max Goplerud, Shiro Kuriwaki, Jens Wiederspohn, Adam Conner-Sax, and Philip Greengard.

Bob Carpenter shared a Stan example to get a flat prior over tables that match specified margins. I think this prior would imply a prior on what the authors call alpha0, the covariance coefficients for the target geography. The method as described in the Modeling Complex Contingency Tables poster seems to focus on point estimates:

In other news, Shiro responded on Twitter to one of my questions about their Application 1 (ACS): “the density plot shown there is an empirical density of 1700+ TREs, where each error is for a non-Southern county.”

I have 2 remaining questions:

  1. What are the covariates w in this ACS example ?
  2. Suppose you also have survey data in Palo Alto. Would this be added to the training tables ?

Their Application 2 asks if lower postratification table reconstruction error improves the downstream MRP:

In the comments last week, Shiro cited related work by Si and Zhou (2021) who propose a method called Bayes-Raking to incorporate known margins into modeling. They found Bayes-Raking was similar to raking in the overall mean but outperformed raking for subgroups:

It would be interesting to directly compare Bayes-Raking to the poster’s SynthMargins !

Survey Statistics: Modeling Complex Contingency Tables

Andrew looped me into an email thread with folks working on poststratification with partial population information. (He knew I’d be interested, see “poststratification without population level information”.)

They pointed to work by Max Goplerud, Shiro Kuriwaki, Jens Wiederspohn, Adam Conner-Sax, and Philip Greengard, which was recently presented as poster at polmeth: Modeling Complex Contingency Tables.

I wish I got to attend, this is cool ! I’ve got questions about Application 1 (ACS):

  1. What are the covariates w?
  2. Suppose you also have survey data in Palo Alto. Would this be added to the training tables ?
  3. How is the table reconstruction error random ? (I see a density plot.)

The poster gives an answer to Andrew’s question in “Mister P when you don’t have the full poststratification table, you only have margins” and “The continuing challenge of poststratification when we don’t have full joint data on the population”:

I’d recommend first imputing a full poststrat table … But then the question is how to do this.

In Application 2, they ask if a better poststratification table reduces error in estimating the outcome using Multilevel Regression and Poststratification (MRP). Presumably this depends on the distribution of the outcome given the poststratification variables. In “toy example for energy balancing weights” I wrote that raking does well “when Y | X1, X2 is additive”.

I’m excited to read the paper and learn more !

Survey Statistics: wanting workflow

Last week Andrew commented that we need a more transparent workflow for survey statistics. So I looked in the new Bayesian Workflow book:

Chapter 19 “Building up to a hierarchical model: Coronavirus testing” is a case study about a 2020 survey that tested n = 3330 residents of Santa Clara County, California for SARS-CoV-2 antibodies (Bendavid et al. 2020a, 2020b). y = 50 people tested positive.

To estimate population prevalence, they want to account for measurement error in the test (“Measurement“) and differences between sample and population (“Representation“), both sides of Groves et al. Figure 2.5:

Image

The measurement error model includes specificity gamma = P[test negative | no disease] and sensitivity delta = P[test positive | disease], which take population prevalence pi to test positivity rate p:

Bendavid et al. (2020a) analyzed their data using gamma = 0.995 and delta = 0.80. But these aren’t known exactly, so Gelman and Carpenter 2020 recommended using priors from previous studies to reflect uncertainty:

They extend these priors to account for multiple studies. They look at sensitivity to these priors (overloaded term here: “sensitivity”). Cool stuff !

In this Survey Statistics series we hadn’t yet talked about models for measurement error in the outcome variable y. We’ve seen measurement error in an adjustment variable X, e.g. recalled vote (“is a mismeasured X better than none at all ?”, “more adventures in mismeasured X”, “more on recalled vote”, “it is (still) the people”).

On the “Representation” side, Bendavid et al. (2020a, 2020b) adjusted for differences between sample and population using weights and encountered the difficulty we saw last week: adjusting for lots of variables can lead to very large weights. Maybe modeling can help (see last week’s “structured MRP to smooth survey weights”).

So Gelman and Carpenter 2020 proposed replacing (19.2) above with

y_i ~ bernoulli(p_i)
p_i = (1-gamma) (1 - pi_i) + delta pi_i

and using Multilevel Regression and Poststratification (MRP):

I thought this story sounded familiar and remembered Andrew talked about it at his Birthday workshop‘s closing talk (this part on YouTube, minute 2 to minute 8).

Andrew said that when his friend tried the method from Gelman and Carpenter 2020 it “crashed and burned” ! There is a new paper about this: Kuh, Kennedy, Chen, and Gelman 2026. They present “a statistical workflow for diagnosing unexpected results in complex models”. These authors have researched a lot about evaluating MRP models (see “should MRP workflow include LOCO-CV ?”). We will have to revisit their new paper later in this series.

Survey Statistics: structured MRP to smooth survey weights

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):

 

Survey Statistics: equivalent models, equivalent weights (locally)

Last month we saw that the Times/Siena Poll is now using energy balancing weights (Huling & Mak, 2024). In a toy example, we saw under which outcome models these weighting methods might do well. I was inspired by Little 2004, who saw under which outcome model the inverse-probability-weighted estimator (a.k.a. Horvitz-Thompson) does well. This explains the HT estimator’s poor performance in Basu’s (1971) elephants example (which I’ve used in my post Basu’s Bears).

From Little 2004:

Andrew’s 2007 “Struggles” paper and the 2026 MrPlew paper go the other way: start with an outcome model and back out the weights. These folks all worked pretty hard. I wondered if I could just use Thomas Lumley’s survey package to get the weights. See my posts “struggles with equivalent weights” and my continued struggles. But in the simulation I used a linear outcome model. It would have been more interesting with a logistic outcome model !

The linear case was handled by Andrew’s 2007 paper. The 2026 MrPlew paper extends this by noticing that the linear equivalent weights are derivatives, and using this to define locally equivalent weights, green highlighting by me:

Then the 2026 MrPlew folks use these equivalent weights to compare the target population to the weighted sample, as an MrP diagnostic. (Unlike the Huling & Mak paper on energy balancing weights, the MrPlew folks look at one covariate function at a time, not the entire covariate distribution at once.) One of their examples is from Lax & Phillips 2009 Gay Rights in the States, which I blogged about for pride. In this example, the implied covariate (im)balance for MrP looks not great ! Here’s their Figure 5b:

As Andrew blogged: “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.” Yes, but here it’s worse than the raw uncorrected balance ! The MrPlew authors caution here that the locally linear approximation may be extrapolating poorly, and that we shouldn’t over-rely on this model check. They also say that “logistic regression generally balances the variance-weighted covariates, but not the covariates themselves”. I am glad they included this illustrative example and I want to think more about it.

Thoughts ?

Survey Statistics: poststratification without population level information

Poststratification uses population data on X to estimate E(Y) via E(E(Y | X, R = 1)), where R = 1 are survey respondents who provide Y and X. When the inner expectation “E” is estimated via Multilevel Regression, this is called MRP. The outer “E” needs p(X), population data on X.

Sometimes we have to estimate the population distributions. We’ve seen a few examples:

  1. “2 flavors of calibration”: Say we have p(X), but we also need p(Z | X), the population distribution of another variable Z. We can estimate p(Z | X, R = 1) using survey data, but nonresponse could make this unreliable. Say we have population data on aggregates p(Z) (e.g. from census tables), then we can logit-shift to anchor to the population aggregate. See Kuriwaki et al. 2024.
  2. “MRPW”: Say we have p(X), but we also need p(W | X), the population distribution of the survey weights W. We can estimate p(W | X, R = 1) using survey data. Then because we assume survey weights are proportional to inverse probability of response 1/p(R = 1 | W, X), we can get what we need by Bayes Rule.
  3. “weights and MRP for voters”: Say we have p(X), but here we need p(X | V = 1), the distribution of X among the population of voters. By Bayes Rule we can get this via p(X) and p(V = 1 | X). The latter can be estimated from population turnout history and vote intent among survey takers.

In all cases, we have some anchor to the population, e.g. via aggregate totals, survey weights, or turnout history. This brings me to Andrew’s post asking 2016 pollsters to poststratify on party ID. We don’t have population aggregates to logit-shift to. But Andrew comments:

party ID is changing much more slowly than the distribution of vote preference, which itself is changing much more slowly than differential nonresponse.

So in our sample party ID (Z) is changing over time quickly, but mostly due to differential nonresponse:

p(Z | t, R = 1) = p(R = 1| Z, t)/p(R =1 | t) * p(Z | t) = differential nonresponse at t * party ID at t

Andrew cites his coauthored paper Reilly et al. 2001, which fits a model to smooth the poststratifying variable Z over time. This won’t help with the component of differential nonresponse that is constant or slowly changing over time, but it prevents the polls from jumping around with every swing in differential nonresponse.

Andrew also cites his coauthored paper The Mythical Swing Voter, which uses 2008 exit poll data on Z to adjust polls from 2012. This also relies on the assumption that Z changes slowly over time, our anchor in the absence of population data.

We’ve been talking about adjusting for party ID. Another variable we might want to adjust for is interest in politics. In “adjusting for interest in politics” we cite Andrew and Gustavo‘s Challenges in Adjusting a Survey That Overrepresents People Interested in Politics. They see an increase over time in interest in politics, which they say could be nonresponse bias (people more interested in politics taking surveys) and/or a change in the population due to increased political polarization.

Survey Statistics: quantifying uncertainty in ranked choice voting polls

We’ve talked about uncertainty in polls (see Margin of Error, Total Margin of Error, Total Margin of Error II) and we’ve talked about ranked data (see exploded logit !). A new paper, Rosenman & Liang 2026, looks at uncertainty in ranked choice voting (RCV) polls.

Recall the multinomial logit model that Train (2009) Chapter 7 calls the exploded logit:

P[ranking Other then Left then Right] = exp(f_Other) / sum_c’ exp(f_c’)   *   exp(f_Left) / (exp(f_Left) + exp(f_Right))

Without covariates, it has only 3 parameters: f_Other, f_Left, f_Right. It makes the independence from irrelevant alternatives (IIA) assumption to go from these 3 parameters to rank probabilities.

In contrast, the multinomial model in Rosenman & Liang 2026 does not make the IIA assumption and has 14 parameters, one for each of 15 possible rankings minus one so they sum to 1:

P[ranking Other then Left then Right] = pi_{Other, Left, Right}

Rosenman & Liang 2026 note that in RCV the election outcome is not expressable as one parameter. Instead, the winner is determined by instant runoff:

  1. If a candidate wins >50% of first choice votes, they win.
  2. Otherwise, the candidate with the least first choice votes is eliminated, and each ballot counts for its top remaining choice. Return to step 1.

Say you use polling data to estimate rank probabilities pi_j for each ranking j. These estimates differ from the true probabilities due to many sources of error (see our favorite Figure 2.5 from Groves et al. shown in quantity vs quality and is a mismeasured X better than none at all ?). Rosenman & Liang 2026 focus on sampling error.

How can we propagate uncertainty about the rank probabilities pi_j to uncertainty about the RCV winner ? If you have draws from the posterior of pi_j, you can do instant runoff on each to get a winner for that draw. This gives win probabilities according to your model and data.

To see the importance of uncertainty in RCV, let’s look at their 2022 Alaska House special election example. With 3 candidates, RCV is determined by 5 margins (see their Lemma 1). Most of these margins are well-identified by the data, but 2 were quite close: Palin vs Begich first choice margin and Peltola vs Palin pairwise margin. They plot these 2 margins in the right panel of Figure 1. The true outcome is the black dot, with sampling uncertainty shown as ellipses around it. For small sample sizes (the biggest ellipse), we see that a plurality of the mass falls into green, where point estimates would declare that Begich wins. Uncertainty quantification would help put this in context, giving all candidates win probabilities around 20-40%, showing the race is difficult to call with such small data.

For details, see Rosenman & Liang 2026.

 

Survey Statistics: toy example for energy balancing weights

Last week we talked about The Big Changes Coming to the Times/Siena Poll:

  1. New weighting variable: support score = E(2024 vote | other X variables).
  2. New weighting method: energy balancing (Huling & Mak, 2024)

Ben Schneider helpfully blogged about energy balancing as well:

Raking and similar calibration methods are based on balancing means or totals for specific variables…The energy balancing method does something different: it calibrates based on an entire multivariate distribution, as measured by an empirical cumulative distribution function (ECDF).

Jared Huling (of Huling & Mak, 2024) helpfully answered questions in the comments. I’m still puzzling over how energy balancing handles empty cells (unsampled regions of the joint covariate space). I need a toy example.

Consider 2 binary variables, so 4 population cells, with known population shares:

       k=0    k=1    total
j=0    .4     .2     .6
j=1    .2     .2     .4
total  .6     .4

Say the sample is missing folks in cell 11:

       k=0    k=1    total
j=0    .5     .3     .8
j=1    .2     0      .2
total  .7     .3

Consider 4 methods:

1. Classical Poststratification: not defined because of division by 0.

2. Raking: match only the margins. Correct when Y | X1, X2 is additive.

       k=0    k=1    total
j=0    .2     .4     .6
j=1    .4     0      .4
total  .6     .4

3. Energy balancing: minimize the Energy-Distance(F_w, F_pop) between the weighted sample distribution of X1, X2 and the population distribution. Correct when Y | X1, X2 is such that nearby cells have similar means.

Say X1 = young/old, X2 = man/woman, Y = percent Democrats, and no old women are sampled.

Raking is correct when additivity holds: old women = young women + (old men − young men)

Energy balancing is correct approximately when: old women = (old men + young women)/2 ?

library(WeightIt)

pop  <- data.frame(X1 = rep(c(0, 0, 1, 1), c(40, 20, 20, 20)),
                   X2 = rep(c(0, 1, 0, 1), c(40, 20, 20, 20)))

samp <- data.frame(X1 = rep(c(0, 0, 1), c(50, 30, 20)),
                   X2 = rep(c(0, 1, 0), c(50, 30, 20)))

dat <- rbind(cbind(pop,  A = 1),
             cbind(samp, A = 0))

W <- weightit(A ~ X1 + X2, data = dat, method = "energy",
              estimand = "ATT", focal = "1",
              dist.mat = as.matrix(dist(dat[, c("X1", "X2")])))

w <- W$weights[dat$A == 0]
tapply(w, interaction(samp$X1, samp$X2), sum) / sum(w)
       k=0    k=1    total
j=0    .381   .309   .69
j=1    .309   0      .309
total  .69    .309

4. MRP: fit a model for Y | X1, X2. The interaction term’s posterior equals its prior, propagating uncertainty around additivity.

Am I understanding this correctly ?

Survey Statistics: Big Changes in the Times/Siena Poll

Yesterday Nate Cohn wrote about The Big Changes Coming to the Times/Siena Poll, with
more details in their poll of Maine.

Say we want to estimate average Platner support in Maine’s likely electorate, E(Y). But we only have survey respondents, R = 1.

The NYT uses survey weights to weight respondents, E(YW | R = 1). In contrast, some pollsters use MRP, fitting a Multilevel Regression model for Platner support, then applying it to the population, E(E_model(Y | X, R = 1)).

Nate discusses 2 Big Changes to how they construct the weights W.

(The polar bear has not yet hiked in ME, but he is training for it. This above is in TN.)

Big Change 1: Support score

A few weeks ago we saw the NYT started weighting on “synthetic 2024 vote”, which is recalled 2024 vote that is validated with the voter file and imputed if needed.

Now they’re also weighting on support score = E(2024 vote | other X variables). Nate explains the motivation:

While a poll can’t weight on dozens of variables, the support score lets us pile a lot of information into a single measure.

This reminded me of the causal inference context, where D’Amour and Franks (2021) “see especially strong performance for propensity weights computed with respect to the prognostic score”, where the prognostic score is E(Y | X, control). In our survey context, this would be a model for Platner support Y. Instead, the NYT use 2024 vote, perhaps for applicability across multiple outcomes Y ?

Big Change 2: Energy balancing

Beyond adding new weighting variables, they’re also changing how they calculate the weights. Nate notes the challenge of weighting on many variables and interactions with typical sample sizes. So they are turning to the R package WeightIt, which implements the energy balancing method from Huling & Mak (2024):

This article introduces a new weighting method, called energy balancing, which instead aims to balance weighted covariate distributions. By directly targeting distributional imbalance, the proposed weighting strategy can be flexibly utilized in a wide variety of causal analyses without the need for careful model or moment specification.

The energy balancing weights do not use outcome Y, but the paper notes that estimates can be improved with a model for Y.

How do energy balancing weights handle the challenge of jointly weighting on many variables with typical sample sizes “without the need for model specification” ?

Survey Statistics: perfect collinearity in the sample but not in the population

In 2019, Andrew blogged about collinearity in Bayesian models. In the comments, he pointed to an example from Bayesian Data Analysis, 2nd edition (BDA2). I think it is a useful example to keep in mind when extrapolating from sample to population. Since folks (like me) may only have BDA3 on their shelf, I thought I’d talk thru it.

Amazon.com: Bayesian Data Analysis, Second Edition (Chapman & Hall/CRC Texts in Statistical Science): 9781584883883: Andrew Gelman, John B. Carlin, Hal S. Stern, Donald B. Rubin: Books

Pretend it is 1980 and we are at the US Census Bureau. We just revamped the occupational coding system, and it’s so much better ! We want 1980-style codes on all our old data that only had 1970-style codes. Let’s trade in our peasant blouses for some shoulder pads.

Say we have double-coded training data (n = 10,000) with:

  • O_1980 = occupation coded in the 1980 coding system
  • O_1970 = occupation coded in the 1970 coding system
  • E = education, either high or low
  • I = income, either high or low

We want to impute O_1980 for the single-coded full dataset (N = 1,000,000) with only O_1970, E, and I.

Consider everyone with the a specific occupation according to the 1970 codes, e.g. Accountants. Say there are 200 accountants in the double-coded training data and they have either high income and high education or low income and low education. They have either OCCUP1 or OCCUP2 according to the 1980 codes.

From BDA2 Table 9.1:

Say we use standard regression software to fit p(O_1980 | O_1970 = Accountants, E, I). It will flag the predictors E and I as perfectly collinear, because in the double-coded training sample, education and income are perfectly correlated.

Suppose you drop education and use only income. The single-coded data actually has some low education and high income folks. The model only uses income, so 90% of them get OCCUP1. But suppose I drop income and use only education. My model only uses education, so only 10% of them get OCCUP1. Who is correct ?

As the authors say:

the truth is that we have essentially no evidence on the split for these units… the occupational split for the ‘E=low, I=high’ units should vary between, say, 90/10 and 10/90. … If some variable should or could be in the model on substantive grounds, then it should be included even if it is not ‘statistically significant’ and even if there is no information in the data to estimate it using traditional methods.

 

Survey Statistics: using MRP in later analyses (pride edition)

Happy pride !

One way I celebrated was by reading Lax & Phillips 2009, Gay Rights in the States: Public Opinion and Policy Responsiveness. It’s on-theme, an example in the MrPlew paper (which I also still need to digest), and I wanted examples of using MRP in later analyses.

Lax & Phillips 2009 studied the relationship between state-level public opinion and state adoption of policies affecting gays and lesbians. Andrew blogged about this work in Nov 2008, Jan 2009, and June 2009 when he wrote:

Fancy statistical analysis can indeed lead to better understanding. Jeff Lax and Justin Phillips used the method of multilevel regression and poststratification (“Mister P”…

The paper’s appendix includes a NYT article and an almost-rainbow-colored plot:

Lax & Phillips 2009 used MRP to estimate state-level public opinion E(y | s). Let

  • y_i = 1 if person i supports laws to protect against discrimination in job opportunities (for example), = 0 otherwise
  • s[i] = state where person i lives, e.g. NY
  • L_s = 1 if state s has laws to protect against discrimination in job opportunities, = 0 otherwise

Their Multilevel Regression (“MR” of MRP) model had race, gender, age, education, state, and poll effects:

They modeled the state effect with state-level predictors (% religious conservatives, % Democratic voters in 2004):

Then they Poststratified (“P” of MRP) to the population:

Then they used the MRP estimate of public opinion as a predictor of whether the state adopts the policy:
Pr(L_s = 1) = logit^-1(a + b * y_s^pred)

From their Figure 1:

Questions:

  1. (How) did Lax & Phillips 2009 incorporate uncertainty in the MRP estimate of public opinion y_s^pred in their later analysis of its effect on policy adoption L_s ?
    Footnote 7 says they incorporated uncertainty for non-MRP estimates:

    if we use an opinion index based on disaggregation instead of MRP estimates, correcting for reliability using an error-in-variables approach (eivreg in Stata)…

  2. Are results sensitive to whether policy adoption L_s is a state-level predictor in the MRP model ?