The statistical methods behind the workshop
Foundation to advanced – the methods at work in the case
This is the second module of the companion curriculum. It works through the statistical methods that the EU-frames case draws on, from the simplest linear model to Bayesian and generalised extensions, using the workshop’s own data throughout: the country-year panel you fit models on, the raw respondent-level records behind it, the specification grid behind the class chart, and the five Multi100 analysts whose results you meet on the day. It follows on naturally from getting started with R and assumes only what that module covers – reading a CSV, a dplyr chain, and the idea of a model formula.
Unlike a textbook chapter, every model on this page is genuinely fitted and every output shown is real. The outputs come in two kinds, and the difference is worth knowing from the start. Models on the committed workshop data – the 270-row country-year panel and the small analyst and grid files – execute when this page is built, so what you see is what the code produces, and you can reproduce every one of them yourself from the shipped files. (The one exception is the Bayesian model of section 8: its data ship with the workshop too, but Stan compiles and samples for a minute or two, so its saved output is displayed instead, with the exact call that made it.) Models on the respondent-level microdata cannot ship: the raw Eurobarometer records are GESIS-licensed and may not be redistributed. Those models were fitted once, on the facilitator’s machine, by the script companion/_model-outputs/fit_models.R in the source repository for the site, and their printed summaries are displayed here exactly as R returned them. The code above each such output is the code that produced it.
1. The predictor lives at the country-year
The original author (OA) modelled individual survey respondents directly. Sixteen Eurobarometer waves pool to 416,698 people interviewed across 27 EU member states between 2004 and 2013, each asked what the European Union means to them. The workshop’s core exercise does not work with those individuals. It works with EUframes_cy.csv, which collapses them to 270 country-year rows – 27 countries observed in each of 10 years – with every framing scale averaged within the cell. For country \(c\) in year \(t\), the cell value is the mean of its respondents’ scores:
\[ \bar{y}_{ct} = \frac{1}{n_{ct}} \sum_{i \,\in\, (c,t)} y_{ict} \]
where \(y_{ict}\) is respondent \(i\)’s framing scale and \(n_{ct}\) the number of respondents in that country-year. Unemployment and GDP growth, which are national figures to begin with, attach directly to each cell.
The data therefore live at three levels – respondents, country-years, and countries – and this module works through all three. Working at the middle level gains two things. The first is practical, and it is the reason that the workshop can exist at all: the raw microdata are GESIS-licensed and cannot be redistributed, whereas the derived country-year averages can be released under an open licence and shipped with the materials. The second is that the country-year is the level at which the predictor genuinely varies. Every respondent interviewed in Germany in 2010 shares one German 2010 unemployment rate, so all the information that the data hold about the alignment of unemployment with EU framing is carried by 270 country-level pairings, not by the hundreds of thousands of individuals. Aggregating to that level discards almost no identifying variation about the relationship of interest.
Averaging within a cell removes all the between-individual variation inside it, so nothing about how framing differs across people of different ages, educations or outlooks survives in the panel. The standard errors no longer reflect 416,698 respondents but 270 cells, which is why the fitted models below carry residual degrees of freedom in the low hundreds rather than the hundreds of thousands. And a relationship among country-year averages need not hold among individuals: reading an aggregate association as if it described people is the ecological inference problem, and the country-year panel cannot speak to the individual-level question at all. Whether that matters is a question about what you are trying to measure rather than a technicality, which is the theme that section 6 returns to.
TipTry it
Read data/EUframes_cy.csv and sum its n_cy column (sum(euframes$n_cy)); you should recover the 416,698 respondents that stand behind the 270 rows. Then run nrow(euframes) to confirm the 270. Holding those two numbers side by side – the people who were interviewed and the rows you actually analyse – makes concrete what aggregation does.
2. Start with the plainest model
Every model in this module is a variation on one idea, so it pays to meet the idea in its plainest form. A linear model proposes that the outcome is a straight-line function of a predictor plus noise. For row \(i\) of the panel:
\[ y_i = \beta_0 + \beta_1 x_i + \varepsilon_i \]
Each symbol has a plain meaning. \(y_i\) is the outcome, here a country-year’s mean cosmopolitan framing score. \(x_i\) is the predictor, its unemployment rate. \(\beta_0\), the intercept, is the expected outcome when the predictor is zero. \(\beta_1\), the slope, is the expected change in the outcome for a one-unit rise in the predictor, and it is the number on which every research question in this case depends. \(\varepsilon_i\) is the error: everything else that moves the outcome, assumed to average out to zero. Ordinary least squares (OLS) estimates the two \(\beta\)s by choosing the line that minimises the sum of squared errors,
\[ (\hat\beta_0, \hat\beta_1) = \underset{\beta_0,\,\beta_1}{\arg\min} \sum_i \left( y_i - \beta_0 - \beta_1 x_i \right)^2 \]
which is what lm() (for ‘linear model’) computes. The name is literal, and a picture makes the point:
Nothing in that geometry is tied to a single predictor. Give the model two,
\[ y_i = \beta_0 + \beta_1 x_{1i} + \beta_2 x_{2i} + \varepsilon_i \]
and the fitted object is no longer a line but a plane floating over the two predictors’ grid. The errors are still the vertical gaps between each point and the fitted surface, OLS still minimises their summed squares, and each slope acquires the reading it keeps for the rest of this module: \(\beta_1\) is the expected change in the outcome for a one-unit rise in \(x_1\) with \(x_2\) held fixed – movement along the plane parallel to one axis. That is all that ‘controlling for’ means. The figure below fits exactly this on the panel, with GDP growth as the second predictor – the very variable that the specification menu can add as a co-predictor – and it is a live three-dimensional object: drag it to rotate, scroll to zoom, and watch the red residual drops stay vertical from every angle.
null
1 With a third predictor the plane becomes a hyperplane that no picture can hold, but the algebra in the arg-min never changes, and neither does the reading of a coefficient: each slope is a comparison with everything else in the model held level. The sections that follow stay with the single-slope version until the menu’s co-predictor choice brings growth back in.
The first argument of lm() is a formula, written outcome ~ predictor and read “outcome as explained by predictor”. The panel was read into euframes at the top of the page, this way:
library(readr)
euframes <- read_csv("../data/EUframes_cy.csv") # ../ because this page lives in companion/Fitting the model and reading it back takes two lines:
first_model <- lm(mcosmo ~ unemp, data = euframes)
summary(first_model)
Call:
lm(formula = mcosmo ~ unemp, data = euframes)
Residuals:
Min 1Q Median 3Q Max
-0.208129 -0.036661 0.008259 0.051255 0.150833
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.516683 0.010133 50.989 <2e-16 ***
unemp -0.001098 0.001055 -1.041 0.299
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.07042 on 268 degrees of freedom
Multiple R-squared: 0.004025, Adjusted R-squared: 0.0003084
F-statistic: 1.083 on 1 and 268 DF, p-value: 0.299This output format recurs, with variations, in every model that this module fits, so it is worth reading slowly once. The Call restates the model. The Residuals block summarises the estimated errors \(\hat\varepsilon_i = y_i - \hat y_i\), the gaps between each observed value and the fitted line. The Coefficients table carries the estimates themselves, one row per \(\beta\), four columns each:
- Estimate is \(\hat\beta\) itself. The
unemprow says a one-point rise in the unemployment rate goes with a fall of about 0.0011 in the cosmopolitan-framing score; the(Intercept)row says a hypothetical country-year with zero unemployment would be expected to score about 0.52. - Std. Error is the standard deviation of the estimate across imagined repeated samples: how precisely the data pin the coefficient down. It is the denominator of everything that follows.
- t value is Estimate divided by Std. Error – how many standard errors the estimate sits from zero. It is unit-free, which is why the workshop can put every submitted result on one chart: whatever the scale of the model, its
tmeans the same thing. - Pr(>|t|) is the p-value: the probability of a
tat least this large in magnitude if the true coefficient were zero. Small values are conventionally read as evidence that the association is not just noise.
Below the table, the residual standard error (0.070 here) estimates the spread of \(\varepsilon\), on 268 degrees of freedom (the 270 observations minus the two \(\beta\)s estimated). Degrees of freedom measure how many independent pieces of information remain for estimating uncertainty, and they matter repeatedly later in this module. Multiple R-squared is the share of the outcome’s variance that the line explains – 0.4% here, essentially nothing – and the F-statistic tests all the predictors at once.
So read as a whole: this first model finds no relationship worth reporting. The slope is small, its t of −1.04 is well inside the range that chance produces, and the p-value of 0.30 confirms it. Yet the workshop’s anchor result, on exactly these 270 rows, is a clearly negative and conventionally significant t = −3.80. Nothing about the data changes between here and there – only the model’s structure does, and the rest of this module works through that structure.
The reason to distrust this simple fit sets the agenda for the rest of this module. The standard errors of OLS assume that the 270 errors \(\varepsilon_i\) are independent and identically distributed, and the estimate itself is only unbiased if the error is uncorrelated with the predictor. Both assumptions are strained here. The 270 rows are not 270 unrelated observations. They are 27 countries observed ten times each, and Denmark’s ten rows share everything durable about Denmark. Those shared national features sit in \(\varepsilon\), they are correlated within country, and – since unemployment levels also differ durably between countries – they are correlated with \(x\) too. The next three sections deal with each consequence in turn: first the standard errors, then the estimate itself.
3. Clustering fixes only the uncertainty
Applied to a panel, the simple regression of section 2 has a proper name: pooled OLS, because it pools all country-years into one undifferentiated cloud. Writing the same equation with the two-index notation of the panel makes the pooling visible – one line for all \(c\) and \(t\), with nothing in the model that records that rows share a country:
\[ \bar{y}_{ct} = \beta_0 + \beta_1 x_{ct} + \varepsilon_{ct} \]
R offers several routes to this fit, and the workshop’s materials use three of them: lm() from base R, feols() from the fixest package (Bergé et al. 2026), and plm() from the plm package (Croissant and Millo 2023), the panel-econometrics toolkit that the published Multi100 analysis used. One habit is worth building before any new setting is introduced. Make different functions agree on the model you already understand. If two functions give the same numbers under settings you can name, you know exactly what each is doing; when their defaults later diverge, you will know which setting moved.
library(fixest)
library(plm)
m_lm <- lm(mcosmo ~ unemp, data = euframes)
m_fx <- feols(mcosmo ~ unemp, data = euframes, vcov = "iid")
m_plm <- plm(mcosmo ~ unemp, data = euframes,
index = c("cntry", "year"), model = "pooling")
tibble(
fitted_with = c("lm()", "feols()", "plm()"),
estimate = c(coef(m_lm)["unemp"], coef(m_fx)["unemp"], coef(m_plm)["unemp"]),
std_error = c(sqrt(vcov(m_lm)["unemp", "unemp"]),
sqrt(vcov(m_fx)["unemp", "unemp"]),
sqrt(vcov(m_plm)["unemp", "unemp"]))
)# A tibble: 3 × 3
fitted_with estimate std_error
<chr> <dbl> <dbl>
1 lm() -0.00110 0.00105
2 feols() -0.00110 0.00105
3 plm() -0.00110 0.00105Identical to every digit, because it is one estimator under three names: vcov = "iid" asks fixest for the classical standard errors that lm() and plm(model = "pooling") report by default. That is the like-for-like baseline. Now change one thing only.
The first repair for the correlated-errors problem keeps the estimate and fixes the uncertainty. Cluster-robust standard errors replace the classical variance formula with one that permits the errors inside each cluster – here, each country – to be correlated however they like:
\[ \widehat{\operatorname{Var}}(\hat\beta) \;=\; (X'X)^{-1} \Big( \textstyle\sum_{g=1}^{G} X_g' \hat\varepsilon_g^{\vphantom{\prime}} \hat\varepsilon_g' X_g^{\vphantom{\prime}} \Big) (X'X)^{-1} \]
The notation is less fearsome than it looks. \(X\) is the matrix of predictor values and \(G\) the number of clusters (27 countries). The formula is a sandwich: the outer ‘bread’ terms are the ordinary OLS ingredients, and the middle ‘meat’ sums each country’s residuals in a way that keeps every within-country correlation intact rather than assuming it away. Only the standard error changes; \(\hat\beta\) is untouched.
library(sandwich)
library(lmtest)
# fixest: one argument
m_cl <- feols(mcosmo ~ unemp, data = euframes, cluster = ~cntry)
coeftable(m_cl) Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.516683404 0.025017527 20.6528565 1.183841e-17
unemp -0.001097567 0.002195328 -0.4999559 6.213113e-01
attr(,"vcov_type")
[1] "Clustered (cntry)"# the same repair applied to the lm() fit via the sandwich package
coeftest(m_lm, vcov = vcovCL(m_lm, cluster = ~cntry, type = "HC1"))
t test of coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5166834 0.0250175 20.653 <2e-16 ***
unemp -0.0010976 0.0021953 -0.500 0.6175
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The two routes agree once the settings are matched: the default small-sample adjustment in fixest corresponds to type = "HC1" in sandwich (and to vcovHC(..., type = "sss") in plm); with type = "HC0" the numbers would drift apart in the trailing digits, which is exactly the kind of discrepancy that the like-for-like habit exists to explain. Substantively, clustering doubles the standard error, and the t falls from −1.04 to −0.50, because the 270 rows were never 270 independent pieces of evidence.
What clustering cannot do is rescue the estimate, and here the pooled estimate itself is the problem. Because the model absorbs no between-country differences, its slope blends two different comparisons: how framing moves when a single country’s unemployment rises over time, and how framing differs between high-unemployment and low-unemployment countries in the cross-section. In these data the two point in opposite directions. Section 5 puts numbers on them: about −0.0056 within countries and +0.0041 between. The pooled slope of −0.0011 is their blend, close to zero and, in many nearby specifications, on the wrong side of it.
That is exactly the pattern that participants see on the class chart. Of the 140 specifications in the cosmopolitan universe (outcome = mcosmo), 41 come out with the ‘wrong’ sign, a positive coefficient running against the claim’s negative direction. Almost all of them sit in one estimator:
| Estimator | Wrong-sign specifications (cosmopolitan universe) |
|---|---|
| Pooled OLS (cluster-robust) | 33 |
| Two-way fixed effects | 8 |
| Country fixed effects only | 0 |
| Random effects | 0 |
| Total | 41 |
Thirty-three of the 41 wrong-sign results are pooled models. The sign flip is overwhelmingly a property of throwing away the panel structure rather than of any other choice on the menu. Clustering is the right way to report uncertainty once you have committed to the pooled estimator; it does nothing to the sign, because the sign is a property of the estimate, not of its standard error.
TipTry it
Read data/spec_grid.csv, keep the cosmopolitan universe (filter(outcome == "mcosmo")), and count the wrong-sign specifications by estimator (filter(!supports_claim) |> count(estimator)). You should recover the 33-and-8 split in the table above. The count falls straight out of the committed grid, so you need not take the table on trust.
4. Compare a country to itself
The second repair goes after the estimate itself. The workshop’s anchor model, the one you reproduce in Task 3, is deliberately constrained: cosmopolitan framing regressed on the unemployment rate, with country and year fixed effects and nothing else.
\[ \bar{y}_{ct} = \beta\, x_{ct} + \alpha_c + \gamma_t + \varepsilon_{ct} \]
Here \(\alpha_c\) is a separate intercept for each country and \(\gamma_t\) a separate intercept for each year. The value of the fixed effects is in what they absorb. The country effects \(\alpha_c\) soak up everything time-constant about a country – its language, institutions, size, and its enduring baseline level of EU enthusiasm – so the unemployment slope is identified only from how a country changes against itself over time, never from the fact that, say, Denmark and Greece differ. The year effects \(\gamma_t\) soak up everything common to all countries in a given year, the EU-wide shocks and the shared arc of the crisis among them. Two-way fixed effects, in one line, compare a country to itself over time, net of whatever moved all of Europe together that year.
Mechanically there are two equivalent ways to fit this, and seeing them agree teaches how the estimator works. The literal way estimates every intercept, putting 26 country dummies and 9 year dummies alongside the slope; this is the ‘least squares dummy variable’ route. The efficient way never estimates them at all: the within transformation subtracts each country’s mean and each year’s mean from both sides,
\[ \left( \bar{y}_{ct} - \bar{y}_{c\cdot} - \bar{y}_{\cdot t} + \bar{y} \right) \;=\; \beta \left( x_{ct} - \bar{x}_{c\cdot} - \bar{x}_{\cdot t} + \bar{x} \right) + \tilde\varepsilon_{ct} \]
where a dot in place of an index means ‘averaged over it’. The demeaning wipes the \(\alpha_c\) and \(\gamma_t\) out of the equation, and OLS on the transformed data recovers exactly the same \(\beta\):
m_lsdv <- lm(mcosmo ~ unemp + factor(cntry) + factor(year), data = euframes)
m_within <- plm(mcosmo ~ unemp, data = euframes, index = c("cntry", "year"),
effect = "twoways", model = "within")
m_twfe <- feols(mcosmo ~ unemp | cntry + year, data = euframes, vcov = "iid")
c(lsdv = coef(m_lsdv)["unemp"], plm = coef(m_within)["unemp"],
feols = coef(m_twfe)["unemp"]) lsdv.unemp plm.unemp feols.unemp
-0.003472985 -0.003472985 -0.003472985 Three functions, one estimator, one number. The fixest summary is the most informative to read:
summary(m_twfe)OLS estimation, Dep. Var.: mcosmo
Observations: 270
Fixed-effects: cntry: 27, year: 10
Standard-errors: IID
Estimate Std. Error t value Pr(>|t|)
unemp -0.003473 0.000913 -3.80372 0.00018214 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
RMSE: 0.034875 Adj. R2: 0.715885
Within R2: 0.058465The slope is −0.00347: within a country, a one-point rise in unemployment goes with a fall of about 0.0035 in the mean cosmopolitan-framing score, net of common yearly shocks. Its t of −3.80 (p = 0.0002) is the workshop’s anchor. The two R-squared lines say different things. The overall adjusted R-squared of 0.72 mostly measures how much the fixed effects themselves explain – country identities are hugely predictive of framing levels – while the within R-squared of 0.058 is the slope’s own contribution: unemployment movements explain about 6% of the within-country variation. Both are true; only the second is about the research question.
The value on record for analyst C6HJR is t = −3.804 (df = 233, N = 270), which is what this model returns and what the workshop reproduces in Task 3.
Two footnotes to the anchor, both of which recur on the specification menu. First, grand-mean centring the predictor shifts the intercept but leaves the slope and its t untouched once fixed effects are present: a centred predictor and its raw form give an identical coefficient, which is exactly why the specification grid keeps only the raw row rather than carrying a redundant centred twin. Second, the standard error above is the classical one, matching the published analysis; ask fixest to cluster by country instead and the same slope’s t becomes −2.31 (p = 0.029). Still negative, still conventionally significant, but a reminder that even after the estimator is settled, the uncertainty convention is a separate, consequential choice. When you compare outputs across analysts, check which convention each is reporting before comparing the ts.
What fixed effects cannot do is the other half of the lesson, and it is the reason that the workshop builds a causal graph before choosing a specification. They are silent on time-varying confounding: anything that changes within a country over the period and moves both its unemployment and its EU framing – a national political scandal, a country-specific austerity programme, a swing in the governing party’s Europe stance – is not absorbed by either set of effects and can bias \(\beta\). They are equally silent on feedback. If EU framing itself feeds back onto economic outcomes, or if both respond to some third time-varying driver, the model has no way to know. A fixed-effects design is a strong answer to time-constant confounders and common shocks, and no answer at all to confounders that move.
5. Within and between are different numbers
Between the pooled and fixed-effects extremes sits a third estimator. Random effects keeps a country term in the model but, rather than estimating 27 separate intercepts, treats the country departures as draws from a distribution:
\[ \bar{y}_{ct} = \beta_0 + \beta_1 x_{ct} + u_c + \varepsilon_{ct}, \qquad u_c \sim N(0, \sigma_u^2) \]
Instead of 27 numbers, the model estimates one: \(\sigma_u^2\), the variance of the country effects. Estimation is by generalised least squares, and its mechanics locate the estimator on the pooled-to-fixed spectrum precisely. Where fixed effects subtract each country’s entire mean, random effects subtract a fraction \(\theta\) of it:
\[ \bar{y}_{ct} - \theta \bar{y}_{c\cdot} = \beta_0 (1 - \theta) + \beta_1 \left( x_{ct} - \theta \bar{x}_{c\cdot} \right) + \text{error}, \qquad \theta = 1 - \sqrt{ \frac{\sigma_\varepsilon^2}{\sigma_\varepsilon^2 + T \sigma_u^2} } \]
with \(T\) the number of years per country. At \(\theta = 0\) this is pooled OLS; at \(\theta = 1\) it is fixed effects; the data place it in between according to how much of the total variance sits between countries.
m_re <- plm(mcosmo ~ unemp, data = euframes,
index = c("cntry", "year"), model = "random")
summary(m_re)Oneway (individual) effect Random Effect Model
(Swamy-Arora's transformation)
Call:
plm(formula = mcosmo ~ unemp, data = euframes, model = "random",
index = c("cntry", "year"))
Balanced Panel: n = 27, T = 10, N = 270
Effects:
var std.dev share
idiosyncratic 0.001612 0.040151 0.337
individual 0.003170 0.056299 0.663
theta: 0.78
Residuals:
Min. 1st Qu. Median 3rd Qu. Max.
-0.1258 -0.0287 0.0014 0.0295 0.1275
Coefficients:
Estimate Std. Error z-value Pr(>|z|)
(Intercept) 0.55283984 0.01324871 41.7278 < 2.2e-16 ***
unemp -0.00525044 0.00081212 -6.4651 1.012e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Total Sum of Squares: 0.50821
Residual Sum of Squares: 0.43965
R-Squared: 0.13492
Adj. R-Squared: 0.13169
Chisq: 41.7974 on 1 DF, p-value: 1.0124e-10The Effects block is the new output to learn. The individual variance (0.0032) is \(\hat\sigma_u^2\), the between-country component; the idiosyncratic variance (0.0016) is \(\hat\sigma_\varepsilon^2\), the within noise; the share column says 66% of the error variance is between-country. The reported \(\theta\) of 0.78 places this fit 78% of the way from pooled towards fixed effects, and the slope duly lands near the country-fixed-effects answer: −0.0053 (z = −6.5) against the −0.0056 that country fixed effects give on their own. (Note that the default random-effects model in plm has country effects only, so its natural comparator is one-way, not two-way, fixed effects.)
Random effects is more efficient than fixed effects when its central assumption holds: that the country effects \(u_c\) are uncorrelated with the predictor. Here that assumption is a strong one, because a country’s unemployment level is very likely tied up with the same durable national features that shape its baseline EU framing – exactly the correlation that fixed effects were designed to sweep away. The Hausman test formalises the comparison. If the assumption held, fixed and random effects would estimate the same quantity and differ only by noise, so a large gap between them is evidence against it. Comparing the two one-way (country) estimators:
fe_country <- plm(mcosmo ~ unemp, data = euframes, index = c("cntry", "year"),
effect = "individual", model = "within")
re_country <- plm(mcosmo ~ unemp, data = euframes, index = c("cntry", "year"),
effect = "individual", model = "random")
phtest(fe_country, re_country)
Hausman Test
data: mcosmo ~ unemp
chisq = 9.8431, df = 1, p-value = 0.001705
alternative hypothesis: one model is inconsistentThe test rejects (p = 0.002): the estimators disagree by more than chance allows, the random-effects assumption is not defensible here, and the fixed-effects estimate is the safer of the two.
There is, though, a more informative response than picking a winner, and it dissolves the fixed-versus-random dilemma rather than settling it. Split the predictor into its two kinds of variation – each country’s average unemployment across the decade, and each year’s departure from that average – and give each its own slope. This is the within–between (or Mundlak) formulation:
\[ \bar{y}_{ct} = \beta_0 + \beta_W \left( x_{ct} - \bar{x}_{c\cdot} \right) + \beta_B\, \bar{x}_{c\cdot} + u_c + \varepsilon_{ct} \]
\(\beta_W\) is the within-country slope, \(\beta_B\) the between-country one. The fit below uses lmer() from lme4 (Bates et al. 2015), the mixed-model workhorse that section 6 introduces properly; the random-intercept term (1 | cntry) plays the role of \(u_c\):
library(lmerTest) # lmer() from lme4, plus degrees of freedom for its t-tests
euframes_wb <- euframes |>
group_by(cntry) |>
mutate(unemp_between = mean(unemp),
unemp_within = unemp - unemp_between) |>
ungroup()
m_wb <- lmer(mcosmo ~ unemp_within + unemp_between + (1 | cntry), data = euframes_wb)
summary(m_wb)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: mcosmo ~ unemp_within + unemp_between + (1 | cntry)
Data: euframes_wb
REML criterion at convergence: -862.4
Scaled residuals:
Min 1Q Median 3Q Max
-3.3294 -0.5742 0.0433 0.6396 3.0258
Random effects:
Groups Name Variance Std.Dev.
cntry (Intercept) 0.003170 0.05630
Residual 0.001612 0.04015
Number of obs: 270, groups: cntry, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.711e-01 3.664e-02 2.500e+01 12.857 1.61e-12 ***
unemp_within -5.645e-03 8.218e-04 2.420e+02 -6.869 5.42e-11 ***
unemp_between 4.143e-03 4.010e-03 2.500e+01 1.033 0.311
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) unmp_w
unemp_withn 0.000
unemp_betwn -0.953 0.000Read the two slopes against everything fitted so far. The within slope, −0.0056, is exactly the country-fixed-effects estimate (fit feols(mcosmo ~ unemp | cntry, ...) and compare), because separating out the between variation is precisely what fixed effects do. The between slope is +0.0041, positive and imprecise, and it tells its own story. Countries with chronically higher unemployment do not frame the EU less cosmopolitanly, and if anything the pattern runs slightly the other way. Now the pooled puzzle of section 3 has an arithmetic answer: the pooled slope of −0.0011 was a weighted blend of −0.0056 and +0.0041, two real patterns with opposite signs, averaged into meaninglessness. The Hausman rejection above is the same fact in test form, firing precisely because \(\beta_W \neq \beta_B\). The estimand lesson arrives early here too: ‘the effect of unemployment on EU framing’ was never one number; there is a within-country number and a between-country number, and a model is a decision about which one you are asking for. (For the fuller argument that this formulation should usually replace the fixed-versus-random contest, see Bell and Jones (2015).)
One practical note completes the estimator sequence. The random-effects estimator is also the reason that the grid holds 840 specifications rather than the full 960. Hold the outcome family at Gaussian and the claim-carrying predictor at unemployment, and the menu’s six remaining axes – six outcomes, two predictor forms, two co-predictor sets, four estimators, five samples and two weighting choices – multiply to 960. Frequency weights do not compose cleanly with the variance-component machinery that a random-effects GLS estimates, so the 120 weighted random-effects cells are omitted as not well defined. That arithmetic is all computable from the committed spec_grid.csv, as are the shares quoted from the grid: 68.7% of specifications in the claim’s direction and 63.9% significant.
6. The error structure decides the evidence
Everything so far has lived at the country-year level. OA’s own study did not. Its anchor result – b = −0.0034, t = −4.03, p < .001 (Table 3, Model 1) – comes from a three-level random-intercept model fitted to the individual respondents. This section moves down to that level and fits the same family of models on the actual microdata: 416,698 respondents, rebuilt from the sixteen GESIS waves exactly as the codebook describes.
NoteHow the outputs in this section were produced
The respondent-level file cannot ship with the workshop (GESIS licence), so every model fitted to it – the multilevel fits below, the weighted fits of section 7, and the generalised models of section 9 – is displayed rather than executed on the page. Each output was produced by the script companion/_model-outputs/fit_models.R in the repository for the site, run on the facilitator’s machine against the rebuilt person-level data. The text shown is R’s output verbatim. (The Bayesian summaries of section 8 come from the same script for a different reason: their model uses the committed panel, but compiling and sampling it takes a minute or two rather than the instant that a live chunk needs.) These are single, deliberately simple fits chosen to teach each method’s output. Nothing beyond them is computed on the microdata anywhere in this workshop. To run them yourself you need the raw waves under your own (free) GESIS registration – the repositories module walks through the registration.
A multilevel model (also called a hierarchical or mixed model) fits the respondents directly while representing the nesting explicitly: respondents sit inside country-years, which sit inside countries. Written at its three levels, with \(i\) indexing respondents:
\[ y_{i(ct)} = \pi_{ct} + \varepsilon_{ict} \qquad\qquad \text{(respondents around their cell mean)} \] \[ \pi_{ct} = \alpha_c + \beta x_{ct} + r_{ct} \qquad \text{(cell means moved by unemployment)} \] \[ \alpha_c = \gamma_0 + v_c \qquad\qquad\quad \text{(countries around the grand mean)} \]
or, substituted into one line:
\[ y_{ict} = \gamma_0 + \beta x_{ct} + v_c + r_{ct} + \varepsilon_{ict} \]
with \(v_c \sim N(0, \sigma_v^2)\), \(r_{ct} \sim N(0, \sigma_r^2)\) and \(\varepsilon_{ict} \sim N(0, \sigma_\varepsilon^2)\). Three variance components now share the noise: \(\sigma_v^2\) for durable country differences, \(\sigma_r^2\) for country-year departures, and \(\sigma_\varepsilon^2\) for person-to-person differences inside a cell. The fixed part (\(\gamma_0 + \beta x_{ct}\)) is shared by everyone; the random part is a sum of draws at each level. This is the random-effects idea of section 5 generalised, and (1 | group) in lmer() notation declares each random intercept. Drawn rather than written, the structure looks like this:
Start deliberately wrong, with only the country level declared:
library(lme4)
person <- readRDS("person_level.rds") # 416,698 respondents; GESIS-licensed rebuild
m_two <- lmer(cosmo ~ unemp + (1 | cntry), data = person)
summary(m_two)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cosmo ~ unemp + (1 | cntry)
Data: person
REML criterion at convergence: 353436.4
Scaled residuals:
Min 1Q Median 3Q Max
-1.65509 -0.93396 0.00253 1.03162 1.80580
Random effects:
Groups Name Variance Std.Dev.
cntry (Intercept) 0.003941 0.06278
Residual 0.136677 0.36970
Number of obs: 416698, groups: cntry, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.474e-01 1.222e-02 2.705e+01 44.80 <2e-16 ***
unemp -4.790e-03 1.953e-04 3.912e+05 -24.52 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
unemp -0.140Read the output in two blocks. The Random effects block reports the estimated variance components (\(\hat\sigma_v^2 = 0.0039\) for countries, \(\hat\sigma_\varepsilon^2 = 0.1367\) for respondents), and the Fixed effects table reads like an lm() coefficient table with one addition – a df column, computed by Satterthwaite’s approximation, saying how many independent pieces of information stand behind each t. And that column exposes the model’s flaw. A slope of −0.0048 with t = −24.5 on 391,200 degrees of freedom looks like overwhelming evidence, but the model has been allowed to treat every respondent as fresh information about unemployment, when every respondent in a country-year shares one unemployment value. The predictor only has 270 distinct values; a model whose df runs to six figures is counting people as if they were countries.
Declare the country-year level and the same data answer differently:
m_three <- lmer(cosmo ~ unemp + (1 | cntry) + (1 | cntry:year), data = person)
summary(m_three)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cosmo ~ unemp + (1 | cntry) + (1 | cntry:year)
Data: person
REML criterion at convergence: 350414.6
Scaled residuals:
Min 1Q Median 3Q Max
-1.77238 -0.96186 0.01493 0.99635 1.89818
Random effects:
Groups Name Variance Std.Dev.
cntry:year (Intercept) 0.001485 0.03853
cntry (Intercept) 0.003759 0.06131
Residual 0.135469 0.36806
Number of obs: 416698, groups: cntry:year, 270; cntry, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.526e-01 1.393e-02 4.441e+01 39.682 < 2e-16 ***
unemp -5.237e-03 8.021e-04 2.531e+02 -6.528 3.61e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
unemp -0.502The slope barely moves (−0.0052), but everything around it does. The standard error quadruples, the t falls from −24.5 to −6.5, and the Satterthwaite df collapses from 391,200 to 253 – the scale of the 270 cells, not of the 416,698 respondents, because the cells are where the genuine predictor contrasts live. The new variance component (\(\hat\sigma_r^2 = 0.0015\)) is small in absolute terms, yet acknowledging it changes the inference fourfold. Compare this fit with the aggregate models of section 5: the three-level person model’s slope of −0.0052 (t = −6.5) is, to two decimals, the aggregate random-intercept estimate. Once the error structure matches the real structure of the data, modelling 416,698 people and modelling their 270 cell means give the same answer about a country-year-level predictor – which is why the workshop can teach the case from the shipped panel.
The genuine additions of the multilevel framework come when you let more than the intercept vary. A random slope lets each country have its own unemployment coefficient, drawn from a distribution whose spread is estimated:
\[ y_{ict} = \gamma_0 + (\beta + b_c)\, x_{ct} + v_c + r_{ct} + \varepsilon_{ict}, \qquad b_c \sim N(0, \sigma_b^2) \]
m_slopes <- lmer(cosmo ~ unemp + (1 + unemp | cntry) + (1 | cntry:year), data = person)
summary(m_slopes)Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: cosmo ~ unemp + (1 + unemp | cntry) + (1 | cntry:year)
Data: person
REML criterion at convergence: 350388.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.76748 -0.95971 0.01417 0.99887 1.89847
Random effects:
Groups Name Variance Std.Dev. Corr
cntry:year (Intercept) 1.201e-03 0.034660
cntry (Intercept) 4.133e-03 0.064291
unemp 3.713e-05 0.006093 -0.44
Residual 1.355e-01 0.368061
Number of obs: 416698, groups: cntry:year, 270; cntry, 27
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.553728 0.015552 20.196050 35.606 < 2e-16 ***
unemp -0.005712 0.001684 18.815257 -3.392 0.00309 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr)
unemp -0.658In the Random effects block, the slope’s standard deviation is 0.0061, larger than the mean slope of −0.0057 itself. Countries do not share one unemployment effect: the model estimates that the effect ranges from clearly negative in some countries to positive in others, with the average sitting where the earlier fits put it (the Corr of −0.44 adds that higher-intercept countries tend to have more negative slopes). And the fixed slope’s df is now about 19 – once each country contributes mainly to its own slope, the effective evidence for the average slope is the 27 countries, not the 270 cells. The df column keeps strict count of where the information genuinely comes from, and it will decide the effect-size story in section 10.
Two of the five Multi100 analysts worked at this level, and they bracket the recorded range. Analyst KEVF1 (in R) reported t = −5.90 with a Satterthwaite-style df of about 136, converting to a partial correlation of −0.452. Analyst PRL47 (in Stata) reported z = −3.59 with the respondent count of nearly 395,000 standing in for df, converting to −0.006. Same data, same model family, same negative direction – and a 75-fold gap in converted effect size produced almost entirely by the df convention you have just watched move from 391,200 to 253 to 19 across three model declarations. The unit of analysis, and the error structure declared around it, is therefore an estimand choice: a decision about what ‘the effect’ refers to and which draws of the world would count as replications, not a software preference. The aggregated models answer a question about country-years. The person-level models appear to answer one about people, but with a predictor that varies only by country-year they mostly re-ask the same question with more elaborate error accounting. Reasonable analysts answering “which question?” differently get genuinely different numbers, before any estimation begins.
Carry the grid arithmetic of section 5 down a level, finally, and what the framework means for the multiverse becomes concrete. Moving to the respondents re-estimates the same 840 cells with more elaborate machinery, and it also widens the menu itself. Person-level covariates that simply do not exist in the aggregates – age, gender, education – become available as a new adjustment fork. The estimator axis changes membership: the multilevel designs of this section stand alongside their fixed-effects and pooled counterparts. The design weights that every respondent carries become usable at last (section 7), and with them one block of cells returns to the menu. The aggregate grid had to drop its 120 weighted random-effects cells because frequency weights do not compose with GLS, but genuine design weights do compose with multilevel models through the pseudo-maximum-likelihood route of the next section, so the weighted multilevel cells are well defined at the person level. Multiply the menu out – six outcomes, two predictor forms, two macro co-predictor sets, four estimators, five samples, two weight choices, and the new demographic-adjustment fork – and the person-level menu holds 1,920 Gaussian cells. Add the outcome-family choice that section 9 introduces, one more binary fork, and the universe doubles to 3,840 person-level specifications. Pooled with the family universe of 1,680 aggregated cells (the 840 of section 5 plus their family twins of section 9), they make a 5,520-specification joint universe in which the level of aggregation is itself a further fork, arguably the most consequential one on the whole menu. None of it has been computed; the arithmetic is there to show that a menu widens geometrically the moment at which the data allow new forks.
7. Weights repair representation
Every Eurobarometer respondent arrives with a number attached: a survey weight. This section is about what those numbers are, how estimation uses them, and why they interact awkwardly with the multilevel framework – a corner of practice where software choices genuinely change results.
Weights exist because samples are not miniature populations. Some people are harder to reach, some decline, some designs deliberately over-sample small groups; the result is that a 60-year-old rural woman may stand for more people than an urban student who was easier to find. A design weight \(w_i\) is, at its cleanest, the inverse of person \(i\)’s probability of ending up in the sample, \(w_i = 1 / \pi_i\). Post-stratification adjustments then align the weighted sample with known population margins. Estimation uses the weights by letting each observation count \(w_i\) times. The weighted mean is the simplest case,
\[ \hat{\bar{Y}} = \frac{\sum_i w_i\, y_i}{\sum_i w_i} \]
and weighted least squares generalises it, choosing coefficients to minimise \(\sum_i w_i \left( y_i - \beta_0 - \beta_1 x_i \right)^2\) so that under-represented kinds of people pull the line with their population share rather than their sample share.
The weighting scheme of the Eurobarometer itself, documented in the weighting overview of the GESIS Eurobarometer Data Service (GESIS Eurobarometer Data Service, n.d.), is a complete worked example of that logic, and the data file for every wave labels its weights plainly. The basic one is W1, “WEIGHT RESULT FROM TARGET” – the w1 of the rebuilt person file. It is a pure post-stratification (in the vocabulary of the Eurobarometer, redressment) weight: an iterative marginal-weighting procedure aligns each national sample with the known margins of its population on sex, age, region (NUTS 2) and size of locality, reproducing each country’s real case count. W1 makes each national sample representative of its own country, and it deliberately does nothing to rescale countries against one another. The companion weights divide the remaining labour. W3 (“WEIGHT GERMANY”) and W4 (“WEIGHT UNITED KINGDOM”) combine the separately drawn East/West German and Great Britain/Northern Ireland samples into united national totals – the same sample splits that the DE and GB recodes in the workshop codebook fold together. The long W5–W30 series (“WEIGHT EU6” through “WEIGHT EU27”, euro zone in and out) add the population-size scaling needed when similarly sized national samples must stand for populations as different as Malta’s and Germany’s, and WEX (“WEIGHT EXTRA POPULATION 15+”) extrapolates results to the actual population aged fifteen and over. The habit that transfers to any survey: find the weighting documentation, identify which weight matches the population that your question is about, and say in the methods which one you used. For a question about the within-country alignment of framing and unemployment, W1 is the natural choice – and it is the natural weights fork for any person-level universe built on this menu.
In R, design-based estimation is the job of the survey package (Lumley 2004). You declare the design once – weights, and where the data supply them, clusters and strata – and then every estimator respects it:
library(survey)
des <- svydesign(ids = ~1, weights = ~w1, data = person) # weights only; no PSU info here
m_svy <- svyglm(cosmo ~ unemp, design = des)
summary(m_svy)
Call:
svyglm(formula = cosmo ~ unemp, design = des)
Survey design:
svydesign(ids = ~1, weights = ~w1, data = person)
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5055600 0.0014713 343.61 <2e-16 ***
unemp 0.0001262 0.0001538 0.82 0.412
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for gaussian family taken to be 0.1398685)
Number of Fisher Scoring iterations: 2This output rewards two readings. As a weighting demonstration: the W1-weighted pooled slope is +0.0001 with t = 0.8 – nothing, and on the wrong side of zero. Weights repair representation, and they have done their job; what they cannot repair is a mis-structured model, and this is the pooled estimator of section 3 again, blending within and between variation across 27 countries. As an output-reading exercise: note that svyglm() reports design-robust (linearised) standard errors as a matter of course – there is no classical-SE option to forget, which is one reason that survey statisticians like the package.
The awkward corner is combining weights with the multilevel framework. The difficulty is conceptual before it is computational. A multilevel model needs to know the selection probability at each level – the country’s probability of inclusion, then the respondent’s probability given the country – but surveys typically publish one all-in weight per person, and there is no unique way to split it. The standard answer is pseudo-maximum-likelihood, weighting each level’s contribution to the likelihood (Rabe-Hesketh and Skrondal 2006), and it comes with a genuine fork: level-1 weights can be scaled in more than one defensible way, the choice can move variance estimates noticeably, and the recommendation in the literature is to report sensitivity across scalings rather than trust any one (Carle 2009).
Because implementations differ in those choices, the same weighted multilevel model can return different numbers in different software, and knowing which convention you are running is part of reporting the model. Three practical facts anchor the R side. First, the weights argument of lmer() is not a survey-weights argument – it treats weights as precision weights, a different quantity entirely, and silently gives design-inconsistent answers if you feed it w1. Second, Stata’s mixed accepts probability weights ([pw = ...]) and implements pseudo-maximum-likelihood with its own scaling options (pwscale()); many of the weighted multilevel results in the applied literature were produced under its conventions. Third, the R package that speaks Stata’s dialect is WeMix (Bailey et al. 2023): it takes one weight column per level and fits by pseudo-maximum-likelihood. In our own comparative survey work it has reproduced Stata’s weighted mixed fixed-effect estimates to three decimals, which makes it the tool of choice when a published Stata result needs re-implementing in R.
Two of the five Multi100 analysts worked in Stata, and their public code turns these abstractions into something you can point at. Analyst PRL47 – the person-level multilevel analyst of section 6 – builds the level-2 cell identifier and fits three-level random intercepts with mixed; these lines are from their Multi100 upload:
egen cquarter=concat(country quarter) // the country-by-quarter cell id
mixed cosmo || country: || cquarter: // three-level random interceptsThe weight is missing from that second line: PRL47’s script carefully assembles a survey weight for every wave (gen weight=w3a4a), then uses it only in descriptive checks (sum cosmo cosmo2 [w=weight]) – it never reaches the model, because carrying it there would mean committing to [pw = ...] and a defensible pwscale() choice, exactly the pseudo-maximum-likelihood fork above. Building a weight is easy; defending its use inside a multilevel model is a decision, and a careful analyst stopping short of that decision is itself instructive. And no estimation options are missing by accident either. A bare mixed estimates by maximum likelihood with model-based standard errors, whereas R’s lmer() defaults to REML, so a line-for-line port disagrees with the original for no substantive reason unless one side’s default is overridden. The correspondence, made explicit for Stata users:
* Stata – PRL47's estimator (maximum likelihood by Stata's default)
mixed cosmo || country: || cquarter:
* what a design-weighted version would have needed
mixed cosmo [pw = weight] || country: || cquarter:, pwscale(size)# the R twin – lmer() defaults to REML, so match Stata's ML explicitly
lmer(cosmo ~ 1 + (1 | country) + (1 | cquarter), data = d, REML = FALSE)
# the design-weighted twin – WeMix, pseudo-maximum-likelihood, one weight per level
mix(cosmo ~ 1 + (1 | country), data = d, weights = c("weight", "w_country"))Back to the data of this case. Supplying a constant 1 at the country level reproduces the common ‘weights at level 1 only’ scenario:
library(WeMix)
person <- person |>
mutate(w_lvl2 = 1) # unit weight at the country level – weighting at level 1 only
m_wemix <- mix(cosmo ~ unemp + (1 | cntry), data = person,
weights = c("w1", "w_lvl2"))
summary(m_wemix)Call:
mix(formula = cosmo ~ unemp + (1 | cntry), data = person, weights = c("w1",
"w_lvl2"))
Variance terms:
Level Group Name Variance Std. Error Std.Dev.
2 cntry (Intercept) 0.003741 0.001394 0.06116
1 Residual 0.136155 0.004246 0.36899
Groups:
Level Group n size mean wgt sum wgt
2 cntry 27 1 27
1 Obs 416677 1 416698
Fixed Effects:
Estimate Std. Error t value
(Intercept) 0.555062 0.015043 36.898
unemp -0.004915 0.001488 -3.302
lnl= -175909.61
Intraclass Correlation= 0.02674 Set this beside the unweighted two-level fit of section 6 and read the comparison carefully, because it contains one more like-for-like trap. The weighted slope, −0.0049, sits next to the unweighted −0.0048: the W1 weights barely move the estimate, which is itself a finding (the framing–unemployment alignment is not an artefact of who is over- or under-represented within national samples). The standard errors, though, are not comparable at a glance – 0.0015 here against 0.0002 there. Most of that gap is not the weights but the variance estimator, since WeMix reports design-linearised (robust) standard errors while lmer() reports model-based ones. Two settings changed between the outputs, and only one of them is the weighting. The habit from section 3 applies at every level of sophistication: before attributing a difference to the interesting setting, equalise the boring ones.
The workshop’s own weighting axis, finally, is a humbler object than any of this. The specification menu offers n_cy weighting – each country-year weighted by the number of respondents behind its average, on the reasoning that a cell mean built from 3,000 people is measured more precisely than one built from 300. That is a precision argument about aggregated cells, not a design weight, and the genuine design weights cannot be reconstructed from the shipped aggregates at all; they live with the microdata, as above. The second Stata analyst in this case, KQXUE, works at exactly this end of the spectrum, and their pipeline shows the whole weighting question dissolving in two lines:
collapse cosmo util, by(cntry year) // unweighted country-year means
xtreg cosmo unemployment, fe // = the within estimator of section 4The collapse averages every respondent equally, so the design weights are gone before any model is fitted, and each cell then enters xtreg with equal say whether 300 people or 3,000 stand behind it – precisely the axis that the menu’s n_cy option toggles. (The Stata–R correspondence continues here too: xtreg, fe after xtset cntry_num year is the within estimator in plm.) The menu’s weighting fork, for the record, moves almost nothing: it contributes less than one per cent of the specification variance in the grid.
8. The same model, stated as a posterior
The multilevel framework has a natural next step. Partial pooling – every random effect shrunk towards the grand mean in proportion to its group’s data – is already half the Bayesian idea, and hierarchical models are the setting where the full version pays for itself most clearly. The shift is one of interpretation before machinery. Frequentist estimation asks “what parameter values make this data likeliest?” and attaches uncertainty via imagined repeated sampling. Bayesian estimation treats the parameters themselves as uncertain quantities, starts from a prior distribution \(p(\theta)\) expressing what is plausible before the data, and updates it through the likelihood into a posterior:
\[ p(\theta \mid y) \;\propto\; p(y \mid \theta)\; p(\theta) \]
Everything you report – point estimates, intervals, probabilities of direction – is then read directly off the posterior distribution. For hierarchical models this buys three concrete things: group-level variances carry their full uncertainty even with few groups (our ten years are few), every derived quantity inherits that uncertainty rather than a plug-in point estimate, and the assumptions that frequentist multilevel models make implicitly become visible, stated priors you can inspect and vary.
Spelled out, the three ingredients are these. The likelihood \(p(y \mid \theta)\) is the model itself – exactly the multilevel equations of section 6, read as a statement of how probable the observed data are under candidate parameter values. The prior \(p(\theta)\) says what is plausible for each parameter before the data arrive; here, mild commitments like “the between-country standard deviation is positive and probably not enormous”. Their product, renormalised, is the posterior \(p(\theta \mid y)\): a full probability distribution over all the parameters jointly, of which any coefficient table is only a summary. For real models that distribution has no closed form, so it is explored by sampling. Markov chain Monte Carlo draws thousands of parameter values whose long-run frequencies match the posterior, and the modern samplers (Stan’s Hamiltonian Monte Carlo, underneath everything below) have made this routine for exactly the hierarchical models that this page has been building. Inference then is the summarising of draws: a posterior mean, an interval containing 95% of them, the share of them below zero. Two books carry a reader from this page to fluency. Gelman and Hill’s classic (2007) is the bridge from single-level regression to multilevel and Bayesian thinking – its statistical chapters have aged well, though its software examples predate Stan and should not be typed in. McElreath’s Statistical Rethinking (2020) is the modern course-style introduction, and its examples exist in a free online translation into brms and tidyverse code (Kurz 2023) that pairs naturally with everything on this page.
The brms package (Bürkner 2017) makes the transition easy: it accepts the same formula syntax as lmer(), translates the model to Stan, and estimates it by Markov chain Monte Carlo (MCMC) – four independent ‘chains’ wandering the posterior, whose agreement is itself the convergence diagnostic. Here is the aggregate random-intercept model of section 5, now with a year intercept alongside, in its Bayesian form. This one uses the committed panel, so you can run it yourself; it compiles for about a minute and samples in seconds:
\[ \bar{y}_{ct} \sim N\!\left( \gamma_0 + \beta x_{ct} + v_c + w_t,\; \sigma^2 \right), \qquad v_c \sim N(0, \sigma_v^2), \quad w_t \sim N(0, \sigma_w^2) \]
library(brms)
m_bayes <- brm(mcosmo ~ unemp + (1 | cntry) + (1 | year), data = euframes,
backend = "cmdstanr", chains = 4, iter = 4000, seed = 2026)
summary(m_bayes) Family: gaussian
Links: mu = identity
Formula: mcosmo ~ unemp + (1 | cntry) + (1 | year)
Data: euframes (Number of observations: 270)
Draws: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
total post-warmup draws = 8000
Multilevel Hyperparameters:
~cntry (Number of levels: 27)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.0627 0.0098 0.0471 0.0847 1.0028 1254 1945
~year (Number of levels: 10)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.0184 0.0063 0.0094 0.0333 1.0020 1862 3555
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.5383 0.0157 0.5077 0.5691 1.0011 1019 1805
unemp -0.0036 0.0009 -0.0054 -0.0018 0.9998 5858 6558
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma 0.0378 0.0018 0.0345 0.0414 1.0005 4855 5142
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).The output’s vocabulary maps onto what you already know, with meanings upgraded. Estimate is now the posterior mean and Est.Error the posterior standard deviation – the analogues of coefficient and standard error. The l-95% / u-95% CI columns are a credible interval, and unlike the confidence interval it means what everyone wants intervals to mean: given model and data, the parameter lies in this range with 95% probability. The genuinely new columns are the diagnostics. Rhat compares the four chains’ answers. At convergence it is 1.00, and values above about 1.01 mean the sampler has not settled and the numbers should not be used. Bulk_ESS and Tail_ESS are effective sample sizes – how many independent draws the correlated MCMC samples are worth. And the Multilevel Hyperparameters block shows the machinery earning its keep: sd(Intercept) for year is estimated at 0.018 with an interval running from 0.009 to 0.033, an explicit statement of how little ten years can say about a between-year variance, where a frequentist fit would hand you a bare point estimate.
Notice where the slope itself landed: −0.0036, between the country-only random-effects estimate (−0.0053) and the two-way fixed-effects anchor (−0.0035), because the year intercepts absorb the common yearly shocks much as the year fixed effects did. That placement is a property of the model, not of the Bayesian machinery, and the direct check is to fit the same model by REML – its exact frequentist twin, live on the committed panel:
m_freq_twin <- lmer(mcosmo ~ unemp + (1 | cntry) + (1 | year), data = euframes)
summary(m_freq_twin)$coefficients Estimate Std. Error df t value Pr(>|t|)
(Intercept) 0.538805859 0.0147467137 54.01694 36.537351 9.024290e-40
unemp -0.003638517 0.0008628051 251.41246 -4.217079 3.455893e-05Identical to the displayed posterior within Monte Carlo error: −0.0036 with a standard error of 0.0009 against a posterior mean of −0.0036 with an Est.Error of 0.0009. With 270 observations and weak priors, Bayesian and frequentist answers agree, as they should, because the data dominate. The traditions come apart, in exactly the direction that this section opened with, where data are thin: compare the twins’ between-year variance statements, a bare REML point estimate against a full posterior interval.
The priors that this model ran with are inspectable, and inspecting them is the habit to build:
prior_summary(m_bayes) prior class coef group resp dpar nlpar lb ub tag
(flat) b
(flat) b unemp
student_t(3, 0.5, 2.5) Intercept
student_t(3, 0, 2.5) sd 0
student_t(3, 0, 2.5) sd cntry 0
student_t(3, 0, 2.5) sd Intercept cntry 0
student_t(3, 0, 2.5) sd year 0
student_t(3, 0, 2.5) sd Intercept year 0
student_t(3, 0, 2.5) sigma 0
source
default
(vectorized)
default
default
(vectorized)
(vectorized)
(vectorized)
(vectorized)
defaultThe defaults in brms are deliberately mild: flat on the slope, a Student-t centred on the median of the data for the intercept, and half-Student-t distributions keeping the variance components positive but otherwise unconstrained. Priors are stated assumptions, which makes them arguable – and re-fitting under different reasonable priors, exactly as the workshop re-fits under different reasonable specifications, is the Bayesian wing of the same robustness discipline. The person-level extension is a one-line change (the section 6 formula, data = person, and patience while 416,698 rows sample), and it leads into the generalised models next.
9. The outcome has a floor and a ceiling
One assumption has gone unexamined since section 2: that a Gaussian linear model is the right shape for this outcome. Strictly, it never was. Every framing scale in this case is bounded between 0 and 1; a Gaussian model has unbounded support, so a straight line extended far enough always predicts impossible values, and nothing inside lm() knows that the walls exist. The assumption is wrong at both levels of the data; what changes between them is how much its wrongness costs, and one pair of pictures makes the argument:
At the aggregate level (right panel) the cell means bunch in the middle of the scale, several residual standard deviations clear of either bound. Fitted values from the models of sections 2 to 5 never approach the walls, the probability that the Gaussian error places beyond them is negligible, and the approximation is safe in practice – safe because of where this particular outcome happens to sit, not because the model respects the bounds. That is a judgement to make consciously, and to re-make whenever an outcome drifts towards its floor or ceiling. At the person level (left panel) the judgement reverses emphatically. A respondent’s score is the share of their mentioned items that are cosmopolitan. A quarter of all respondents sit at exactly 0 and another quarter at exactly 1, and a Gaussian model treats those pile-ups as ordinary noise around a middle that describes almost nobody. The generalised linear model (GLM) family fixes the shape problem by modelling a transformed mean:
\[ g\!\left( \mathbb{E}[y_i] \right) = \eta_i = \beta_0 + \beta_1 x_i \]
The linear predictor \(\eta_i\) stays a familiar straight line; the link function \(g\) maps the outcome’s constrained mean onto the unconstrained line, and a distribution suited to the outcome’s actual support replaces the Gaussian.
The classic member is logistic regression, for binary outcomes, with the log-odds link. Collapsing our outcome to ‘mentions at least one cosmopolitan item’:
\[ \log \frac{p_i}{1 - p_i} = \beta_0 + \beta_1 x_i, \qquad p_i = \Pr(y_i = 1) \]
m_logit <- glm(I(cosmo > 0) ~ unemp, family = binomial(), data = person)
summary(m_logit)
Call:
glm(formula = I(cosmo > 0) ~ unemp, family = binomial(), data = person)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 1.262274 0.008360 151.00 <2e-16 ***
unemp -0.026735 0.000837 -31.94 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 481477 on 416697 degrees of freedom
Residual deviance: 480478 on 416696 degrees of freedom
AIC: 480482
Number of Fisher Scoring iterations: 4Reading a GLM’s coefficients takes one extra step, because the estimate lives on the link scale. The slope of −0.027 is a change in log-odds per unemployment point; exponentiating gives an odds ratio of \(e^{-0.027} \approx 0.974\), and the quick ‘divide by four’ rule bounds the probability change at about −0.027 / 4 ≈ −0.7 percentage points per unemployment point at most. The z value replaces t (the dispersion is fixed, not estimated) and deviance replaces R-squared as the fit currency. And one number should now trigger the section 3 alarm: z = −31.9 is the pooled-independence illusion again, 416,698 respondents counted as if each were fresh evidence about 270 unemployment values. The generalised and multilevel machineries are two orthogonal upgrades, and this outcome needs both at once.
For a bounded continuous outcome the classical tool is beta regression, which models a variable on the open interval (0, 1) with a logit-linked mean and a precision parameter \(\phi\). Its defect for our variable is those endpoint masses: the beta distribution has no probability at exactly 0 or 1, and the traditional dodge of nudging the data inward is exactly the kind of silent choice that this workshop exists to make visible. Ordered beta regression (Kubinec 2023) repairs it at the root. One model carries ordered cutpoints deciding whether an observation is exactly 0, in between, or exactly 1, with a beta distribution for the middle, and one coefficient vector governs all three parts, so the slope keeps a single interpretation.
Both of the day’s inferential traditions implement it. The frequentist route is glmmTMB (Brooks et al. 2017), a fast maximum-likelihood engine for GLMs with random effects – here with the full three-level structure of section 6, on all 416,698 respondents:
library(glmmTMB)
m_ordfreq <- glmmTMB(cosmo ~ unemp + (1 | cntry) + (1 | cntry:year),
family = ordbeta(), data = person)
summary(m_ordfreq) Family: ordbeta ( logit )
Formula: cosmo ~ unemp + (1 | cntry) + (1 | cntry:year)
Data: person
AIC BIC logLik -2*log(L) df.resid
683220.8 683297.4 -341603.4 683206.8 416691
Random effects:
Conditional model:
Groups Name Variance Std.Dev.
cntry (Intercept) 0.015782 0.1256
cntry:year (Intercept) 0.006162 0.0785
Number of obs: 416698, groups: cntry, 27; cntry:year, 270
Dispersion parameter for ordbeta family (): 9.28
Conditional model:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 0.183815 0.028629 6.421 1.36e-10 ***
unemp -0.011789 0.001665 -7.080 1.44e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The output assembles pieces you now recognise: variance components for the two grouping levels (on the latent logit scale, so not directly comparable with those of section 6), a dispersion parameter (\(\phi\) = 9.3, the beta precision), and a conditional-model slope of −0.012 on the logit scale with z = −7.1 – the divide-by-four rule translates it to at most −0.3 percentage points of cosmopolitan share per unemployment point, the same order as every Gaussian estimate of it, now from a model whose predictions cannot leave the unit interval. Note the z of −7 rather than −32. The random effects are doing the clustering work that the bare logistic model ignored.
The Bayesian route is ordbetareg, Kubinec’s own implementation of the model on top of brms (Kubinec 2023). The exhibit below runs on a random subsample of 20,000 respondents – the full-data posterior would sample for hours, and the purpose here is the output’s shape, not a production estimate:
library(ordbetareg)
set.seed(2026)
person_sub <- slice_sample(person, n = 20000)
m_ordbayes <- ordbetareg(cosmo ~ unemp + (1 | cntry), data = person_sub,
true_bounds = c(0, 1), backend = "cmdstanr",
chains = 4, iter = 2000, seed = 2026)
summary(m_ordbayes) Family: ord_beta_reg
Links: mu = identity
Formula: cosmo ~ unemp + (1 | cntry)
Data: data (Number of observations: 20000)
Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
total post-warmup draws = 4000
Multilevel Hyperparameters:
~cntry (Number of levels: 27)
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sd(Intercept) 0.1305 0.0209 0.0962 0.1782 1.0051 647 1264
Regression Coefficients:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept 0.1508 0.0316 0.0894 0.2107 1.0163 464 1126
unemp -0.0078 0.0020 -0.0118 -0.0038 0.9998 4370 3540
Further Distributional Parameters:
Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
phi 9.2926 0.1278 9.0450 9.5445 1.0009 3555 2674
cutzero -0.9643 0.0170 -0.9976 -0.9314 1.0000 3123 3302
cutone 0.7709 0.0085 0.7543 0.7878 1.0002 3138 3064
Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).Everything in it combines section 8 with what came before, and it rewards the discipline of reading diagnostics first. The intercept’s Rhat of 1.016 sits above the 1.01 bar, so its chains have not fully settled. In a production analysis the response would be to sample longer before using that parameter. The slope, by contrast, is clean (Rhat at 1.00, effective sample sizes in the thousands), and its posterior – −0.008 on the logit scale, credible interval −0.012 to −0.004 – agrees with the full-data glmmTMB estimate while being less precise on a twentieth of the data. The family’s extra distributional parameters close the loop: cutzero and cutone are the ordered cutpoints governing the exact-0 and exact-1 masses, and phi, the beta precision, posts 9.29 where the maximum-likelihood fit in glmmTMB put it at 9.28 on all 416,698 respondents. The two estimation traditions fit the same model family and return the same answer.
The reason to know this model family goes back to the workshop’s central finding: the fork that moves the results in this case is the outcome’s operationalisation, and how you model the outcome’s distribution is part of that fork. Treat the family as one more binary fork, Gaussian versus beta, both defensible for a bounded scale, and every universe it touches doubles. For the aggregated grid the doubling is computed, not hypothetical. The committed data/spec_grid_family.csv holds the family universe of 1,680 specifications – the canonical 840, imported verbatim, plus a beta twin of each – built by companion/_model-outputs/build_family_grid.R from the committed panel and grid alone. The file sits on the specification menu as a first-class axis, and one label travels with it everywhere: the family axis is a workshop extension, not part of the Multi100 menu.
A second extension works the predictor side of the same lesson. Task 2 of the Multi100 project fixed unemployment as the measure of poor economic performance, but Task 1 shows that three of the five analysts, left free, chose GDP measures instead – and OA’s own operationalisation used both. The committed data/spec_grid_full.csv therefore adds GDP growth as a second, claim-carrying predictor: raw form, alone or beside unemployment as a co-predictor, 840 cells across both outcome families, taking the full universe to 2,520 (the family universe of 1,680 rows plus these 840 growth rows). Only the raw form is fitted, deliberately. Affine transforms – centring, z-standardising – rescale a coefficient without moving its t, p or partial correlation, so on an r-ranked curve they would only ever duplicate a point that is already there. That is exactly why the canonical grid keeps a single raw-and-log pair per specification rather than four affine variants, and it is also part of why the predictor-form fork’s share of specification variance is essentially zero. The growth universe returns its own result: the direction agrees with the claim about as often as it does for unemployment (70.5%), but significance falls below half (46.9%).
The twins are plain logit-link beta regressions on the cell means, well defined because the 270 means sit strictly inside the unit interval, as the right panel above showed, with none of the endpoint masses that force the ordered-beta repair at the person level. betareg carries the pooled and fixed-effects designs, with the fixed effects as country and year dummies and the pooled errors cluster-robust by country. beta_family() from glmmTMB carries the random-effects design, and the weighted random-effects exclusion of section 5 is mirrored, 120 cells per family. A beta model reports a z statistic rather than a t, and that changes nothing downstream. The section 10 conversion \(r = z / \sqrt{z^{2} + \mathit{df}}\) reads a z exactly as the record read analyst PRL47’s, so every twin lands on the common axis unchanged; the grid also evaluates each twin’s significance against the same t reference that its Gaussian sibling uses, a deliberate uniformity that matters only in the third decimal place at 100-plus df. Both halves can be checked straight from the committed file:
read_csv("../data/spec_grid_family.csv", show_col_types = FALSE) |>
group_by(family) |>
summarise(specifications = n(),
claim_direction = mean(supports_claim))# A tibble: 2 × 3
family specifications claim_direction
<chr> <int> <dbl>
1 beta 840 0.710
2 gaussian 840 0.687The twins leave two facts worth carrying away. The Gaussian half posts the canonical shares, 68.7% of specifications in the claim’s direction and 63.9% significant; the beta half posts 71.0% and 63.1%. The beta twin of the workshop’s anchor – mcosmo, raw predictor, no co-predictor, two-way fixed effects, all years, unweighted – estimates a logit-scale slope of −0.014 with z = −4.075, converting to a partial correlation of −0.258 against the Gaussian anchor’s −0.242. Same direction, still significant, slightly stronger. And the wrong-sign anatomy of section 3 reappears line for line. In the cosmopolitan universe each family holds exactly 41 wrong-sign specifications, 33 of them pooled and 8 two-way fixed effects. Run the ‘Try it’ count from section 3 on the family grid and the split repeats twice over. The family fork moves the numbers without changing the lessons.
The person-level doubling, unlike the aggregated one, is arithmetic rather than computation. Carrying the family fork down a level meets one substitution that the person-level data force: the person-level scales place real mass on exact zeros, so the beta likelihood is unavailable there and a logit-link fractional-response binomial would stand in as each cell’s twin (reverse-coding stays exact, since the flip works through the same logit identity). That takes the person-level universe to 3,840 specifications and the joint universe to 5,520, none of it computed. What the arithmetic teaches is that universes grow geometrically, one defensible fork at a time, so the discipline is to state in advance which forks you consider defensible and why, rather than to run the largest universe you can afford – which is what a preregistration is for.
10. Put every result on one axis
The five Multi100 analysts fitted different models, in different software, at different levels of analysis, so their raw coefficients are not comparable: a Stata pooled coefficient, an R fixed-effects coefficient and a person-level multilevel coefficient are not the same kind of number. The discipline with a century of practice at exactly this problem is meta-analysis, and its machinery is worth meeting properly before seeing what Multi100 borrowed from it – not least because it is a model that this page has already fitted.
Meta-analysis starts from \(K\) estimates \(\hat\theta_k\), each with a known sampling variance \(v_k\) (from its standard error). The equal-effects model assumes one true value \(\theta\) behind them all and pools by inverse-variance weighting, so precise studies count for more:
\[ \hat\theta_{\text{EE}} = \frac{\sum_k \hat\theta_k / v_k}{\sum_k 1 / v_k} \]
The random-effects model drops the one-true-value assumption and lets the true effects themselves vary between studies:
\[ \hat\theta_k = \mu + u_k + e_k, \qquad u_k \sim N(0, \tau^2), \quad e_k \sim N(0, v_k) \]
Here \(\mu\) is the average true effect, \(\tau^2\) the between-study variance, and \(e_k\) each study’s sampling error. Read that equation beside section 6 and the correspondence is exact: a random-effects meta-analysis is a two-level multilevel model – studies as level-2 units, random intercepts \(u_k\) – with a single twist, that the level-1 variances \(v_k\) are not estimated but known, imported from each study’s standard error. The two traditions differ in little else. Heterogeneity statistics then quantify the between-study spread: \(\tau\) on the effect’s own scale, Cochran’s \(Q\) as a test, and \(I^2\) as the share of total variability that is real between-study difference rather than sampling noise.
For correlations, the conventional pooling scale is Fisher’s z transformation, \(z = \operatorname{atanh}(r)\), whose sampling variance depends only on sample size, \(v = 1/(n - 3)\). The five analysts’ recorded partial correlations are in the committed analysts5.csv, so the standard tool – the metafor package (Viechtbauer 2010) – can run on them live:
library(metafor)
analysts <- read_csv("../data/analysts5.csv", show_col_types = FALSE) |>
mutate(zi = atanh(r), vi = 1 / (n - 3))
ma_re <- rma(yi = zi, vi = vi, data = analysts, method = "REML", slab = analyst_id)
ma_re
Random-Effects Model (k = 5; tau^2 estimator: REML)
tau^2 (estimated amount of total heterogeneity): 0.0472 (SE = 0.0347)
tau (square root of estimated tau^2 value): 0.2172
I^2 (total heterogeneity / total variability): 99.97%
H^2 (total variability / sampling variability): 3250.05
Test for Heterogeneity:
Q(df = 4) = 31765.7952, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub
-0.2481 0.0991 -2.5035 0.0123 -0.4423 -0.0539 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Run as a demonstration, this output is worth every line; run as an analysis, it would be seriously wrong, and both halves are the lesson. Take the heterogeneity block first, where \(I^2\) reaches 99.97%. In a genuine meta-analysis that number would say the studies barely share a phenomenon; here it says something cleaner, because we know that all five results come from the same data and claim. Sampling error explains essentially none of the spread between the five analysts – the spread is analytic variation, choices about outcome, estimator, level and df, which is precisely the subject of the multiverse showing up in a heterogeneity statistic. The pooled estimate, by contrast, should not be believed even as arithmetic, because meta-analysis assumes independent studies and five reanalyses of one dataset are the textbook violation. The weighting adds its own distortion – the two person-level analysts enter with \(v_k\) built from six-figure respondent counts, so the equal-effects pooled value (run method = "EE" yourself and it returns −0.17) is mostly their voice, inheriting exactly the df problem that section 6 dissected. The Multi100 design takes the point: it borrowed the standardisation step of meta-analysis, but deliberately not its pooling step, reporting the distribution of results rather than an average pretending to independence.
The native picture of meta-analysis, the forest plot, makes the same case at a glance:
forest(ma_re, atransf = tanh, refline = 0,
xlab = "partial correlation r (back-transformed from Fisher's z)",
header = c("Analyst", "r [95% CI]"))
One brief, five defensible readings
The five reanalyses were not free-for-alls. Task 2 of the Multi100 project fixed both ends of the question – use the positive non-materialist (cosmopolitan) form of EU framing as the outcome, use the unemployment rate as the measure of poor economic performance, and disregard the EU-policy and party-politicisation contextual variables of the paper – and left everything between them open. What each analyst did with that freedom, read off their public code:
| Analyst | Data as analysed | Outcome construction | Predictor | Controls | Estimator | SEs and df |
|---|---|---|---|---|---|---|
| 018OL | 420 country-by-wave cell means, EU-27 | mean of the 5 items, rescaled 0–1 | unemployment 15–74, raw | none | crossed random intercepts: country, year, EB wave | Satterthwaite, 396.9 |
| C6HJR | 270 country-years | cosmopolitan share, EU-frames panel | unemployment, centred | none | two-way FE: country + year | conventional, 233 |
| KEVF1 | 209,826 respondents | sum of the 5 items (0–7) | unemployment, raw | none | nested random intercepts: year within country | Satterthwaite, 135.6 |
| KQXUE | 270 country-years | cosmopolitan share | unemployment 20–64, raw | none | one-way country FE | conventional, 242 |
| PRL47 | 394,575 respondents | cosmopolitan share, no-opinion coded 0 | unemployment 25–64, raw + yearly change | gender, age, education, occupation, rural, crisis dummy, log GDP, growth | three-level: respondent in country-quarter in country | Wald z, no df |
That brief, though, was the second act of the two-phase Multi100 design (Aczel et al. 2026). Task 1, months earlier, was free: analysts “could conduct and report as many analyses as they wished, but they had to draw a single conclusion”. What freedom produced belongs beside the standardised table:
| Analyst | Free-phase predictor | Free-phase outcome | Recorded verdict |
|---|---|---|---|
| 018OL | GDP per capita (level, Eurostat) | positive framing, 7-item mean | evidence for the claim |
| KEVF1 | GDP growth (Eurostat) | positive framing, 7-item sum | evidence for the claim |
| PRL47 | log GDP + quarterly growth + crisis dummy | three positive-framing shares | evidence for the claim |
| KQXUE | GDP growth and unemployment (Eurostat) | cosmopolitan + utilitarian shares | evidence for the claim |
| C6HJR | GDP growth and unemployment (World Bank) | positive framing, 7-item share | evidence for the claim |
Three of the five never used the unemployment rate until the instructions named it, and most were not using the non-materialist outcome that the instructions fixed. Yet all five recorded the same categorical Task-1 verdict, because Task 1 asked for a conclusion, not a number – no standardised Task-1 effect sizes exist anywhere in the public Multi100 data. The instructions therefore did not trim covariates from five similar models. They created the common question that the standardised results answer. Read together, the two phases separate two robustness claims that are easily conflated: the five free analyses agreed in words, and the five standardised estimates, finally comparable, disagreed across a 75-fold range. This workshop runs the arc in reverse – participants start from the standardised baseline and open one axis at a time – so the day’s freedom arrives as a menu rather than a blank page.
Every row of the standardised table is a defensible reading of the same brief, and the table explains the numbers that puzzle first-time readers. The Ns differ because the units differ. KQXUE and C6HJR collapsed respondents to 270 country-year means; 018OL kept the survey wave as its own grouping level, giving 420 country-by-wave cells; KEVF1 and PRL47 modelled respondents directly, and even their two Ns disagree for pipeline reasons – KEVF1 merged an unemployment series with full country coverage only from 2009 (2004–2008 exists only for France), so complete-case fitting halved the stacked file to 209,826, while PRL47’s EU-27 filter and respondent-level controls left 394,575. The dfs then differ because each model keeps its own reckoning of where the information about a country-year predictor lives: fixed-effects OLS counts observations minus parameters, the Satterthwaite approximation estimates an effective df (397 on 420 cells, but 136 on 209,826 respondents, because the predictor varies only across country-years and the approximation knows it), and Stata’s mixed reports an asymptotic z with no residual df at all. The conversion machinery below turns those choices into the numbers on the common axis.
The standardisation step is this. Every recorded result is converted to a partial correlation from just two numbers, the test statistic and its degrees of freedom:
\[ r = \frac{t}{\sqrt{t^{2} + \mathit{df}}} \]
Feed it any model’s t (or z) and df, and it returns a number between −1 and 1 that a Stata pooled model, an R fixed-effects model and every class submission can be read on together. This is the exact conversion that the workshop’s report_result() helper applies before a dot lands on the multiverse chart, and it is why Task 2 asked every analyst for a statistic and its degrees of freedom in the first place. The closely related Cohen’s \(d\) follows from the same \(r\):
\[ d = \frac{2r}{\sqrt{1 - r^{2}}} \]
The same pair of numbers also carries the significance verdict, and the conversion is worth stating alongside the effect-size one because the two run in opposite directions. A t statistic’s two-sided p-value is the probability of a result at least as extreme when the true effect is zero – in R, 2 * pt(-abs(t), df). A z statistic runs the identical computation on the standard normal, 2 * pnorm(-abs(z)), which is simply the t-distribution with its df sent to infinity. That is literally where PRL47’s missing df entry went. The df matters only at the margins here: with 20 df the 5% threshold is |t| = 2.09, and by a hundred df it has effectively converged on the normal’s 1.96. What matters far more is what each conversion does with the volume of data. The p-value asks whether the estimate stands out from zero given how much data went in, so at any fixed underlying effect it shrinks as n grows; the partial correlation divides the sample size back out. Feed both conversions the same t and df and they can disagree spectacularly about which results are ‘big’, which is exactly what the five analysts show:
| Analyst | Tool | Model unit | Statistic | df | p (two-sided) | N | Partial r |
|---|---|---|---|---|---|---|---|
| 018OL | R | aggregate | t = −1.21 | 396.9 | .229 | 420 | −0.060 |
| C6HJR | R | country-year (270) | t = −3.80 | 233 | < .001 | 270 | −0.242 |
| KEVF1 | R | person-level multilevel | t = −5.90 | 135.6 | < .001 | 209,826 | −0.452 |
| KQXUE | Stata | country-year (270) | t = −7.17 | 242 | < .001 | 270 | −0.419 |
| PRL47 | Stata | person-level multilevel | z = −3.59 | – | < .001 | 394,575 | −0.006 |
Every analyst found a negative effect, agreeing on direction, and four of the five clear the conventional significance bar decisively. Read the last two columns against each other, though, and the verdicts come apart. PRL47’s z of −3.59 is unambiguously significant (p < .001) while converting to a partial correlation of −0.006 – significance carried by 394,575 respondents, not by the size of the effect. 018OL, the one non-significant row (p = .23), converts to an r ten times PRL47’s. A p-value and a standardised effect are answers to different questions, and neither substitutes for the other. Yet the partial correlations run from −0.006 to −0.452, a factor of about 75, and the two person-level analysts sit at the opposite extremes. You have now seen the whole mechanism from the inside: KEVF1’s conversion divides through a group-scale df of about 136, of the kind that the three-level and random-slope fits of section 6 produce; PRL47’s runs the respondent count through the denominator, of the kind that the two-level fit produced. The same negative relationship, converted under two df conventions, becomes either a substantial correlation or a rounding error. The habit that the workshop insists on follows: read the degrees of freedom before reading the effect size, because a converted number is only as meaningful as the df behind it.
The common axis is also what makes the workshop’s central result legible. Decomposing the 840-specification grid to ask which analytical fork moves the estimate most gives an unambiguous answer:
| Analytical fork | Share of specification variance |
|---|---|
| Outcome (which framing scale) | ≈ 49.7% |
| Estimator + predictor form + co-predictor + weights (combined) | < 2% |
The choice of outcome – which of the four framing dimensions, or their positive and negative composites, you treat as ‘EU framing’ – accounts for roughly half of all the variation across the grid. Everything else that the menu lets you vary, the estimator choices of sections 3 to 6 included, together moves the result by under two per cent. Two qualifications keep that number readable. The decomposition is computed on the grid as estimated, before the claim-alignment convention that the workshop’s charts apply for display. Part of the outcome share is therefore the polarity split between positive and negative framings: real variation in what the specifications estimate, which the aligned chart removes from the picture but not from the decomposition. Recomputed on the claim-aligned scale in the facilitator’s extended analysis of the same grid, the ranking inverts, with the estimator fork rising to about 58% and the outcome fork falling to roughly 11–12%. Once every outcome points the same way, the pooled-versus-within contrast of section 5 becomes the biggest mover. The outcome-family fork of section 9, meanwhile, barely registers: across the family universe of 1,680 rows the Gaussian-versus-beta choice accounts for under 1% of claim-aligned specification variance, and the two family curves all but coincide.
The lesson that participants leave with holds under every recount: deciding what ‘EU framing’ means – which scale, and which polarity – sets what a model estimates before the machinery of sections 3 to 6 gets its say, and modelling the outcome’s distribution matters far less than choosing the outcome itself. For what a specification curve can and cannot conclude once all those results are on one axis, the workshop’s reading on the specification curve picks up exactly here.
Sources and provenance
This module is written fresh for the workshop rather than adapted from an existing lesson, and every number in it is computed, not quoted. The models of sections 2 to 5, the outcome-scale figure and family-grid summary of section 9, and the meta-analysis of section 10 execute on this page from the committed files (EUframes_cy.csv, analysts5.csv, spec_grid_family.csv, cosmo_dist.csv). The country-year construction and item definitions are in the data codebook. The person-level and Bayesian outputs of sections 6 to 9 were produced by companion/_model-outputs/fit_models.R in the source repository for the site, run against the person-level rebuild of the sixteen GESIS waves. The raw microdata are GESIS-licensed and never travel with the workshop – the repositories module explains the registration route. The Stata excerpts in section 7 are from analysts PRL47’s and KQXUE’s public Multi100 materials; the Task-1 free-phase recipes and recorded verdicts in section 10 come from the same analysts’ public code and the public dataset of the Multi100 project (Aczel et al. 2026). The wrong-sign counts, the 840-versus-960 arithmetic, and the claim-direction and significance shares come from spec_grid.csv. The fork-importance decompositions (outcome ≈ 49.7% of specification variance as estimated, estimator ≈ 58% on the claim-aligned scale) come from the facilitator’s extended analysis of the same grid, and the 3,840- and 5,520-specification universes of sections 6 and 9 are menu arithmetic, not computed analyses. The exhibit fits on this page are illustrations of each method; nothing beyond them is computed on the microdata.
The methods and software cited above are gathered in the reference list below: the panel estimators (Bergé et al. 2026; Croissant and Millo 2023), within–between modelling (Bell and Jones 2015), the multilevel workhorse (Bates et al. 2015), survey-weighted estimation (Lumley 2004; Rabe-Hesketh and Skrondal 2006; Carle 2009; Bailey et al. 2023; GESIS Eurobarometer Data Service, n.d.), Bayesian estimation and its textbooks (Bürkner 2017; Gelman and Hill 2007; McElreath 2020; Kurz 2023), the generalised models (Brooks et al. 2017; Kubinec 2023), and meta-analysis (Viechtbauer 2010).
References
Aczel, Balazs, Barnabas Szaszi, Harry T. Clelland, et al. 2026. “Investigating the Analytical Robustness of the Social and Behavioural Sciences.” Nature 652 (8108): 135–42. https://doi.org/10.1038/s41586-025-09844-9.
Bailey, Paul, Blue Webb, Claire Kelley, Trang Nguyen, and Huade Huo. 2023. WeMix: Weighted Mixed-Effects Models Using Multilevel Pseudo Maximum Likelihood Estimation. https://american-institutes-for-research.github.io/WeMix/.
Bates, Douglas, Martin Mächler, Ben Bolker, and Steve Walker. 2015. “Fitting Linear Mixed-Effects Models Using Lme4.” Journal of Statistical Software 67 (1): 1–48. https://doi.org/10.18637/jss.v067.i01.
Bell, Andrew, and Kelvyn Jones. 2015. “Explaining Fixed Effects: Random Effects Modeling of Time-Series Cross-Sectional and Panel Data.” Political Science Research and Methods 3 (1): 133–53. https://doi.org/10.1017/psrm.2014.7.
Bergé, Laurent R., Kyle Butts, and Grant McDermott. 2026. “Fast and User-Friendly Econometrics Estimations: The R Package Fixest.” arXiv Preprint arXiv:2601.21749, ahead of print. https://doi.org/10.48550/arXiv.2601.21749.
Brooks, Mollie E., Kasper Kristensen, Koen J. van Benthem, et al. 2017. “glmmTMB Balances Speed and Flexibility Among Packages for Zero-Inflated Generalized Linear Mixed Modeling.” The R Journal 9 (2): 378–400. https://doi.org/10.32614/RJ-2017-066.
Bürkner, Paul-Christian. 2017. “Brms: An R Package for Bayesian Multilevel Models Using Stan.” Journal of Statistical Software 80 (1): 1–28. https://doi.org/10.18637/jss.v080.i01.
Carle, Adam C. 2009. “Fitting Multilevel Models in Complex Survey Data with Design Weights: Recommendations.” BMC Medical Research Methodology 9 (1): 49. https://doi.org/10.1186/1471-2288-9-49.
Croissant, Yves, and Giovanni Millo. 2023. Panel Data Econometrics in R: The Plm Package. https://cran.r-project.org/web/packages/plm/vignettes/A_plmPackage.html.
Gelman, Andrew, and Jennifer Hill. 2007. Data Analysis Using Regression and Multilevel/Hierarchical Models. Cambridge University Press.
GESIS Eurobarometer Data Service. n.d. Weighting Overview (Standard and Special Eurobarometer). Accessed July 17, 2026. https://www.gesis.org/en/eurobarometer-data-service/data-and-documentation/standard-special-eb/weighting-overview.
Kubinec, Robert. 2023. “Ordered Beta Regression: A Parsimonious, Well-Fitting Model for Continuous Data with Lower and Upper Bounds.” Political Analysis 31 (4): 519–36. https://doi.org/10.1017/pan.2022.20.
Kurz, A. Solomon. 2023. Statistical Rethinking with Brms, Ggplot2, and the Tidyverse: Second Edition. https://bookdown.org/content/4857/.
Lumley, Thomas. 2004. “Analysis of Complex Survey Samples.” Journal of Statistical Software 9 (8): 1–19. https://doi.org/10.18637/jss.v009.i08.
McElreath, Richard. 2020. Statistical Rethinking: A Bayesian Course with Examples in R and Stan. 2nd ed. CRC Texts in Statistical Science. Taylor; Francis, CRC Press.
Rabe-Hesketh, Sophia, and Anders Skrondal. 2006. “Multilevel Modelling of Complex Survey Data.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 169 (4): 805–27. https://doi.org/10.1111/j.1467-985X.2006.00426.x.
Viechtbauer, Wolfgang. 2010. “Conducting Meta-Analyses in R with the Metafor Package.” Journal of Statistical Software 36 (3): 1–48. https://doi.org/10.18637/jss.v036.i03.