Fake data, planted coefficients, and what they teach about real analyses
This is the simulation module of the companion curriculum, and it is optional. Every other quantitative page on this site analyses data that came from the world; this module analyses data that came from a recipe we wrote ourselves, so that for once the answer is known before any model runs. That inversion, known as fake-data simulation, is a standard habit in applied statistics. Regression and Other Stories uses it constantly, on the argument that fitting models to data whose coefficients you planted yourself is how you learn what your models can and cannot recover (Gelman et al. 2020); the same habit anchors classroom practice (Gelman and Vehtari 2024) and the workflow literature, where experimentation on simulated data is one of the named components of model-building (Gelman et al. 2026).
The module teaches from a dataset built for the purpose. EUframes_person_full.csv is a simulated person-level companion to the EU-frames case – 41,660 synthetic respondents across the real 270 country-year cells, with the real macro variables attached. It carries the original analysis rather than a cut-down version of it: all four framing dimensions and both composites, plus every individual variable that the original author (OA) used – sex, age, education as both years and bands, an eight-category occupational position, urbanisation, and left-right self-placement.
First, no Eurobarometer respondent appears in it. The raw microdata behind the case are GESIS-licensed and cannot be redistributed in any form, including as samples or extracts – the repositories module explains that regime – so a person-level dataset that this site can legitimately ship has to be synthetic. What the licensed file supplied, read once and locally, is a set of summaries: marginal distributions, a correlation matrix, and the coefficients of a model of the survey instrument. Those are aggregate statistics of the kind that published appendices carry, and they are committed as _model-outputs/joint_*.csv and _model-outputs/item_*.csv.
Second – and this is the pedagogical point rather than a limitation – the recipe is known. A data-generating process is the explicit rule by which every value was produced. Real data never comes with its rule attached, whereas simulated data is nothing but the rule, so any model fitted to it can be judged against an answer that is already known.
The generator is _model-outputs/simulate_person_joint.R, a committed script worth reading, and it works in two stages that answer two different questions: who is in the sample, and what they say.
Stage one: who. A respondent’s characteristics are one draw from a latent multivariate normal whose correlation matrix was estimated on the real respondents. Each observed variable is a monotone transform of one coordinate of that draw: sex and employment status by a threshold, urbanisation by cutpoints, age, education and left-right through their quantile functions. Any marginal distribution you like can be reproduced this way, and every pairwise association comes with it.
The associations matter as much as the marginals. An earlier version of this twin drew age, sex and education independently, which is the natural thing to do and is wrong: in the real respondents, age correlates with being unemployed at -0.24 and with education at -0.14. A regression coefficient is a partial association, so it depends on how the variables covary; get that wrong and the coefficients come out wrong however carefully each variable’s own effect was set.
Stage two: what they say. In the real survey a respondent hears a list of things that the EU might mean and names those that apply, and their dimension scores are the shares of their mentions falling to each dimension. The generator simulates that instrument rather than the four scores. Each of the thirteen items has a latent propensity
\[ z_{ij} = \alpha_{j(ct)} + \mathbf{x}_i^{\top}\boldsymbol\beta_j + \varepsilon_{ij}, \qquad \text{item } j \text{ mentioned when } z_{ij} > 0, \]
with \(\mathbf{x}_i\) the respondent’s characteristics, \(\alpha_{j(ct)}\) a shift for their country-year, and the thirteen residuals correlated with one another. That is a multivariate probit, and each of its three parts does a different job: \(\boldsymbol\beta\) carries who says what, the correlated residuals carry the fact that people who name one thing tend to name its neighbour, and \(\alpha\) carries the country-year.
Simulating the instrument rather than the scores is what gives the shares their exact zeros and ones, makes them sum to one automatically, and lets every dimension respond on its own terms.
bcoef |> slice_max(abs(estimate), n = 6) |>
mutate(estimate = signif(estimate, 3))# A tibble: 6 × 3
item term estimate
<chr> <chr> <dbl>
1 it6 edu3high 0.397
2 it5 edu3high 0.321
3 it7 edu3high 0.302
4 it9 edu3high -0.264
5 it5 sesstudent 0.244
6 it3 edu3high 0.241Occupation dominates that list, which is why the twin carries it. An earlier version left occupation out of the answer stage entirely – respondents had an occupation, but it had no bearing on what they said – and the result reproduced the macro slopes while getting two thirds of the coefficients in the model wrong. Section 2 is where that was found.
The country-year shift is solved, not set. A dimension score is a share of a respondent’s own mentions, so raising the cosmopolitan items lowers every other dimension even with its own items untouched. There is no formula for the shift that produces a given set of cell means, so for each of the 270 country-years the generator solves four shifts by Newton iteration against a fixed reference sample, until the expected dimension means match those in the panel to within 1e-4.
sim |>
slice_sample(n = 6) |>
select(cntry, year, age, edu3, ses, urban, cosmo, util, comm, lib, unemp) |>
arrange(cntry, year)# A tibble: 6 × 11
cntry year age edu3 ses urban cosmo util comm lib unemp
<chr> <dbl> <dbl> <fct> <fct> <fct> <dbl> <dbl> <dbl> <dbl> <dbl>
1 AT 2008 46 mid manager rural 0 0 0 0 4.13
2 BE 2005 45 low retired large_town 0.5 0.5 0 0 8.44
3 IT 2004 30 mid manager small_town 1 0 0 0 7.87
4 LT 2012 52 mid unemployed rural 0 0 0 0 13.4
5 LU 2010 53 low white_collar large_town 0.6 0.4 0 0 4.36
6 SI 2010 36 mid manual rural 0 0 0 1 7.24# The share columns are stored rounded to six decimals, so four of them can
# drift from their sum: the identities need a tolerance, not ==.
identities <- sim |>
filter(cosmo + util + comm + lib > 0) |>
summarise(
`dimensions sum to 1` = all(abs(cosmo + util + comm + lib - 1) < 5e-6),
`pos + neg = 1` = all(abs(pos + neg - 1) < 5e-6),
`share exactly 0` = mean(cosmo == 0),
`share exactly 1` = mean(cosmo == 1))
identities# A tibble: 1 × 4
`dimensions sum to 1` `pos + neg = 1` `share exactly 0` `share exactly 1`
<lgl> <lgl> <dbl> <dbl>
1 TRUE TRUE 0.197 0.224Among respondents who mention anything, 20% score exactly zero on cosmopolitan framing and 22% score exactly one – the latter because someone who mentions one item and finds it cosmopolitan scores 1 by construction. This is faithful to the real instrument, and it is the concrete reason that beta regression, which needs the outcome strictly inside the unit interval, cannot be used at person level while it can be used on country-year means. The statistical methods module states that constraint; here you can verify it.
The first thing to do with data whose recipe you know is to ask a model to find what the recipe planted. There are two such questions here, and they need keeping apart.
At cell level the planted value is the coefficient from the panel itself, because the generator solves each country-year until its expected dimension means match those in the panel. So a model fitted to the cell means of the twin should land where the same model fitted to the real panel lands.
map(planted$outcome, \(v) {
cf <- coef(lm(anchor(v), data = sim_cells))
tibble(outcome = v, sim_unemp = cf[["unemp"]], sim_growth = cf[["growth"]])
}) |>
list_rbind() |>
left_join(planted, by = "outcome") |>
select(outcome, planted_unemp, sim_unemp, planted_growth, sim_growth) |>
mutate(across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 6 × 5
outcome planted_unemp sim_unemp planted_growth sim_growth
<chr> <dbl> <dbl> <dbl> <dbl>
1 cosmo -0.00356 -0.00412 -0.000242 -0.000899
2 util -0.000678 -0.000186 0.00191 0.00169
3 comm 0.00524 0.00516 -0.00163 -0.00174
4 lib 0.000105 0.000485 0.00000514 -0.0000412
5 pos -0.00424 -0.00431 0.00167 0.000788
6 neg 0.00534 0.00564 -0.00162 -0.00179 Read the substance before the arithmetic. Rising unemployment pushes framing away from both cosmopolitan and utilitarian terms and towards communitarian ones. Growth works differently again, lifting utilitarian framing most – which makes sense, since that dimension is built from the economic-prosperity and social-protection items. That pattern belongs to the real panel, carried into the twin rather than manufactured by the simulation.
Now the arithmetic. The simulated coefficients sit near the planted values but not on them, and the reason lies in the construction: each cell mean in the twin is an average over a tenth of the real respondents, so it is noisier than the corresponding mean in the panel. That is a property of the design – not a defect – and a single dataset could never show it. Only repeated draws from the same recipe can.
ggplot(sweep, aes(sim_unemp, outcome)) +
geom_point(alpha = 0.45, colour = "#00205b", size = 1.6) +
geom_point(data = planted, aes(planted_unemp, outcome),
colour = "#c8102e", shape = 124, size = 8) +
labs(x = "estimated unemployment slope", y = NULL) +
theme_minimal(base_size = 12)
planted_long <- planted |>
select(outcome, sim_unemp = planted_unemp, sim_growth = planted_growth) |>
pivot_longer(-outcome, names_to = "slope", values_to = "planted")
sweep_summary <- sweep |>
pivot_longer(c(sim_unemp, sim_growth), names_to = "slope",
values_to = "estimate") |>
summarise(draws = n(), mean_estimate = mean(estimate),
sd_estimate = sd(estimate), .by = c(outcome, slope)) |>
left_join(planted_long, by = c("outcome", "slope")) |>
mutate(gap_in_ses = (mean_estimate - planted) / (sd_estimate / sqrt(draws)))
sweep_summary |>
arrange(slope, outcome) |>
mutate(across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 12 × 7
outcome slope draws mean_estimate sd_estimate planted gap_in_ses
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
1 comm sim_growth 19 -0.00173 0.000753 -0.00163 -0.555
2 cosmo sim_growth 19 0.000154 0.00103 -0.000242 1.67
3 lib sim_growth 19 -0.000184 0.000484 0.00000514 -1.7
4 neg sim_growth 19 -0.00191 0.000987 -0.00162 -1.26
5 pos sim_growth 19 0.00187 0.00115 0.00167 0.74
6 util sim_growth 19 0.00171 0.000505 0.00191 -1.73
7 comm sim_unemp 19 0.00599 0.000703 0.00524 4.68
8 cosmo sim_unemp 19 -0.0042 0.00101 -0.00356 -2.74
9 lib sim_unemp 19 0.00016 0.000632 0.000105 0.383
10 neg sim_unemp 19 0.00615 0.00101 0.00534 3.48
11 pos sim_unemp 19 -0.00497 0.00115 -0.00424 -2.78
12 util sim_unemp 19 -0.000772 0.000414 -0.000678 -0.983The last column is the test: each remaining gap as a multiple of the standard error of the mean across draws. 8 of the 12 fall inside two standard errors, and the ones outside are not scattered. Every unemployment slope misses in the same direction, with the estimate from the twin larger in magnitude than the planted value by 16 per cent at the median. The median rather than the mean, because the planted slope of libertarian framing is nearly zero and any absolute error there is an enormous proportion of it – an average would report a number driven entirely by the one quantity with almost nothing in it to estimate.
That is a real bias in the generator, and it is shown here because of how it had to be found. No single draw reveals it. Fit the anchor to the committed twin and regress its cell means against those in the panel, and the unemployment coefficient on the difference comes back at t = −0.12 – nothing at all. One draw’s sampling noise is the same size as the bias, so a sweep of this size is what separates it from zero.
Its source is a stated limitation of the construction rather than a mistake in it. The generator gives each country-year four numbers, one shift per dimension, where the real data has thirteen, one rate per item. A dimension-level shift captures as little as 37 per cent of the cell-to-cell variation in some items, and what is lost shows up twice over: as a mentioning rate about 1.6 points below the one in the panel, and as this amplification.
The obvious repair does not work, and anyone tempted to try it should know that first. Giving each cell its own thirteen item rates makes the marginals exact – mentions per respondent land within 0.05 per cent of target – and makes the individual-level validation worse, from 39 coefficients agreeing to 32. In the real data a cell’s item rates reflect its context and the mix of people in it together; the twin draws its people from pooled marginals, so forcing the real rates onto a different mix pushes the difference into the country-year term, and it re-emerges as distorted covariate coefficients. Fixing it properly would mean reproducing each cell’s demographic composition as well, which is a much larger undertaking than the gain justifies.
Where that leaves the twin: unbiased for the individual-level model it was validated against, and biased by about a sixth for cell-level unemployment slopes, with the reason understood and the cost of removing it judged too high.
At individual level the planted values are not those of the panel. They are whatever the real respondents give under the same model, so finding out whether the twin reproduces them means fitting that model on both files. That comparison is what licenses the twin to be used at all, and it lives in _model-outputs/twin_validation.csv:
valid |>
summarise(coefficients = n(),
`within 1 SE` = sum(abs(z) < 1),
`within 2 SE` = sum(abs(z) < 2),
largest = round(max(abs(z)), 2))# A tibble: 1 × 4
coefficients `within 1 SE` `within 2 SE` largest
<int> <int> <int> <dbl>
1 64 46 60 2.8760 of 64 agree within two standard errors of the difference. The individual-level page works through what those coefficients say.
That validation is also how two real failures in this generator were found, and neither was visible any other way. The first version left occupation out of the answer stage, so a model containing occupation and education split the association between them differently from the real data – two thirds of the coefficients wrong, while every marginal, every cell mean and every correlation checked out. The second carried education as a straight line through a relationship that the published model reads in three bands, and returned about half of each band contrast.
Any simulated dataset can pass every check you think to write and still fail the analysis you built it for, and the only check that catches that is running the intended analysis on both files.
Open the generator and change SIM_FRACTION from 0.10 to 0.40. It is the only lever on how many respondents each cell gets, and nothing about the recipe’s coefficients changes – only how precisely a cell mean can be measured.
Do not judge the result on one rerun: a single draw gives one estimate, and one estimate cannot tell you whether the spread has narrowed. Sweep instead – ten seeds at each setting, with $SIM_SEED naming the draw and $SIM_SWEEP collecting the estimates – then compare the standard deviation of each set of ten. Quadrupling the cell sizes should roughly halve it, which is the square-root law arriving as an experimental result rather than a formula. That is a designed simulation study in miniature, and a power analysis is built the same way.
Point $SIM_OUT at a scratch directory unless you mean to overwrite the committed twin, and note that at 0.40 the file is four times the size.
On the grid as estimated, the choice of outcome – which framing scale you treat as ‘EU framing’ – moves results more than any other analytical decision. The twin lets you see part of why, because the dimensions are shares of one instrument and cannot move independently.
planted |>
select(outcome, unemp = planted_unemp, growth = planted_growth) |>
mutate(ratio_to_cosmo = unemp / unemp[outcome == "cosmo"],
across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 6 × 4
outcome unemp growth ratio_to_cosmo
<chr> <dbl> <dbl> <dbl>
1 cosmo -0.00356 -0.000242 1
2 util -0.000678 0.00191 0.191
3 comm 0.00524 -0.00163 -1.47
4 lib 0.000105 0.00000514 -0.0293
5 pos -0.00424 0.00167 1.19
6 neg 0.00534 -0.00162 -1.5 On unemployment the signs fall exactly along the line that the day’s claim-alignment convention describes, both positive framings declining as unemployment rises and both negative ones gaining. Growth is what breaks the convention. There the cosmopolitan dimension moves with the communitarian one and utilitarian framing with the libertarian, so the split runs between whichever dimension is gaining and whichever are giving way, and which is which depends on the predictor. Libertarian framing is the weak case in that reading, since its growth coefficient is under one per cent of the communitarian one in size and what the two share is a sign rather than a magnitude. What claim alignment actually rests on is narrower and does hold exactly here: pos and neg partition the mentions, so neg = 1 − pos for every mentioner and reverse-coding the negative composite is exact arithmetic.
The magnitudes also span a wide range, so an analyst choosing util and one choosing comm are estimating genuinely different quantities from one instrument. That is the structure of the panel itself rather than a simulation artefact, and it is a substantial part of why the outcome fork carries so much specification variance. The multiverse module takes that argument further, because it is precisely where Auspurg’s measurement-nonequivalence objection bites.
A country-year’s unemployment rate could move EU framing in two ways. It might be that unemployed people frame the EU differently, so a high-unemployment cell holds more of them and its mean shifts with the mix – a compositional effect. Or living through high unemployment might change how everyone frames the EU, whatever their own situation – a contextual effect.
The country-year panel cannot separate them, because it contains no individuals. This twin can, because it carries a calibrated occupational structure: being unemployed has its own estimated effect on what a respondent says, and the share of the respondents in a cell who are unemployed rises with the national rate by a measured amount.
Both halves matter, and the second is the one that gets assumed. A large individual effect spread over a group that barely grows moves an average hardly at all, and the arithmetic that puts the two together tells you more than either coefficient alone. The individual-level page does that calculation on this twin, and finds the compositional path accounts for a few per cent of the country-level association.
What belongs here is the methodological point rather than the answer. An earlier version of this twin manufactured the compositional channel: it gave respondents an invented unemployment indicator whose prevalence was a fixed multiple of the national rate, precisely so this section would have something to display. Two things were wrong with that. Nothing about it was measured – the prevalence rule and the effect were both chosen. And it fused two constructs that live at different levels, a labour market’s unemployment rate and a person’s employment status, which quietly asserted a mechanism that the case never claimed and leaked into every macro coefficient that the twin planted.
The general rule matters more than any particular exhibit of it. A simulation can only demonstrate what was put into it. If the mechanism on display is one you wrote down yourself, the demonstration establishes that your code works, not that the world does – which is why simulated data can be used to check whether an estimator recovers an answer, but never to supply one.
Age, sex, education, occupation and urbanisation all carry into the twin the effects they have in the real respondents, which is what lets it support the published model rather than a simplified one. Whether they arrive intact is a question with a checkable answer.
valid |>
filter(outcome == "cosmo",
!term %in% c("(Intercept)", "munemp", "mgrowth")) |>
select(term, real, twin, z) |>
mutate(across(where(is.numeric), \(x) signif(x, 3))) |>
arrange(desc(abs(real)))# A tibble: 13 × 4
term real twin z
<chr> <dbl> <dbl> <dbl>
1 edu3high 0.117 0.102 -2.62
2 sesstudent 0.0916 0.086 -0.718
3 sesmanager 0.0639 0.0599 -0.567
4 edu3mid 0.0584 0.0462 -2.5
5 sesunemployed -0.035 -0.0278 0.984
6 seswhite_collar 0.0303 0.0232 -1.09
7 sesself_employed 0.0238 0.0218 -0.259
8 urbanrural -0.0192 -0.0186 0.111
9 urbansmall_town -0.0117 -0.0131 -0.321
10 sesretired 0.00663 0.0145 1.13
11 seshouse_person -0.00257 0.0121 1.93
12 female 0.00186 0.00932 1.97
13 mage -0.0012 -0.000973 1.47 Education is the largest individual effect in the file: graduates sit about 0.12 higher on cosmopolitan framing than those who left school by fifteen, and the twin returns 0.10. The z column reads each gap against the standard error of the difference, so a value inside two means the twin reproduces that coefficient as well as a second sample of its size would.
Notice how much more precisely these individual coefficients are pinned down than the macro ones in section 2. Forty thousand respondents carry information about education; 270 cells carry information about unemployment. The same dataset is richly informative about one question and thin on the other, which is the design lesson that section 9 makes explicit.
Build your own placebo. Take any person-level variable in the file, shuffle it within country-year, and refit the model above with the shuffled version in place of the real one. Its coefficient should land near zero with an interval covering it, and if you repeat that twenty times about one run in twenty will not, which is the five-per-cent rule doing exactly what it advertises. It is far more instructive to meet a false positive where you are certain there is nothing to find than to read about one.
The italics on person-level are doing real work. unemp, growth and bailout hold one value for a whole country-year, so shuffling them inside a cell puts every value back where it was and hands you the real coefficient with the real p-value, looking for all the world like a placebo that failed. That is not a quirk of the shuffle: a variable that does not vary within the unit you are permuting cannot be permuted there at all, which is the same aliasing that makes a null impossible to build at the level at which a macro predictor lives.
The statistical methods module tells a story it can only display, because it rests on the licensed file: a naive analysis that ignores clustering reports an enormous t on hundreds of thousands of degrees of freedom, while the multilevel model reports a modest one on a couple of hundred – because the predictor varies only at the country-year level. On the twin you can run the contrast.
m3 <- lmer(cosmo ~ unemp + (1 | cntry) + (1 | cntry:year), data = sim)
naive <- lm(cosmo ~ unemp, data = sim)
m3_par <- parameters::model_parameters(m3, ci_method = "satterthwaite")
df_tbl <- tibble(
model = c("naive OLS, clustering ignored", "three-level random intercepts"),
estimate = c(coef(naive)[["unemp"]], fixef(m3)[["unemp"]]),
se = c(sqrt(diag(vcov(naive)))[["unemp"]],
m3_par$SE[m3_par$Parameter == "unemp"]),
t = estimate / se,
df = c(df.residual(naive), m3_par$df_error[m3_par$Parameter == "unemp"]))
# The two rows differ in the estimate as well as in the standard error, so it is
# worth having the within-country slope to hand: it is the one that the
# multilevel fit leans towards, while the naive fit blends it with the
# between-country one.
cells_w <- sim_cells |>
left_join(count(sim, cntry, year), by = c("cntry", "year"))
B_WITHIN <- coef(lm(cosmo ~ unemp + factor(cntry), data = cells_w,
weights = n))[["unemp"]]
df_tbl |> mutate(across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 2 × 5
model estimate se t df
<chr> <dbl> <dbl> <dbl> <dbl>
1 naive OLS, clustering ignored -0.000555 0.000441 -1.26 41700
2 three-level random intercepts -0.00564 0.00108 -5.24 251The naive row claims forty thousand independent witnesses; the multilevel row counts the few hundred effective ones, because respondents inside a country-year share one unemployment value and largely echo their cell. Person-level sample sizes advertise precision that the design does not possess whenever the predictor lives at a higher level. The model that respects the structure pays for that respect in degrees of freedom, and here in the estimate too: the naive fit is a different weighted average of between-country and within-country variation, landing at -5.55^{-4} against the multilevel -0.00564, which sits near the within-country slope of -0.00549. Neither of the two is the planted coefficient, because both leave out the year effects that the anchor carries. This is also why, in Part 1 of the workshop, two Multi100 analysts fitting person-level models to the same claim could reach such different strengths of evidence.
Survey weights are usually taught as a procedure. The interesting question is when they change anything, and simulated data can answer it because the distortion that the weight corrects can be built in deliberately. The sample here was drawn with a fieldwork distortion – younger respondents over-selected – and w1 is the exact reciprocal of that selection factor, normalised to mean 1 within each country-year.
# Cell sizes are drawn in proportion to each cell's real n_cy, so a person-level
# mean over the respondents in the twin estimates the respondent-weighted value
# in the panel. The average of the 270 cell means is a different quantity and
# would make the weight look like it was doing harm.
POP_COSMO <- weighted.mean(panel$mcosmo, panel$n_cy)
ERR_UNW <- mean(sim$cosmo) - POP_COSMO
ERR_WTD <- weighted.mean(sim$cosmo, sim$w1) - POP_COSMO
tibble(
quantity = c("mean age", "mean cosmopolitan share"),
unweighted = c(mean(sim$age), mean(sim$cosmo)),
weighted = c(weighted.mean(sim$age, sim$w1),
weighted.mean(sim$cosmo, sim$w1)),
population = c(AGE_REF, POP_COSMO)
) |>
mutate(across(where(is.numeric), \(x) round(x, 3)))# A tibble: 2 × 4
quantity unweighted weighted population
<chr> <dbl> <dbl> <dbl>
1 mean age 43.1 48.1 48.0
2 mean cosmopolitan share 0.479 0.471 0.5The achieved sample is about four years too young, and weighting recovers the population age almost exactly. Because age genuinely affects framing in the recipe, that composition error propagates into the outcome: unweighted, the cosmopolitan share sits -0.021 above the population value, and weighting removes about -38 per cent of that. It re-balances who is counted, and it matters exactly to the extent that the mis-counted people differ on the thing being measured.
The population column is itself a choice. A person-level mean over the respondents in the twin is aiming at the respondent-weighted value in the panel, not at the average of its 270 cell means, because cells enter the twin in proportion to how many real respondents they held. The two differ by 0.007, which is larger than the whole effect of the weight, so a table benchmarked against the wrong one would show weighting making the estimate worse. Choosing the estimand – the quantity that a statistic is actually estimating – decides which direction the lesson runs in.
The distortion here is constant across cells, so it shifts levels without touching the relationship of interest – the common case, and the reason that analysts often find weighting changes their descriptive statistics more than their coefficients. It stops being the common case the moment the distortion varies with the predictor.
In the generator, set SIM_DISTORT_SLOPE to 0.01. The youth over-selection then grows with a cell’s unemployment rate rather than sitting at one strength everywhere, and because age is correlated with framing in the estimated structure, the distortion now bends the macro slope itself instead of only shifting levels. That is the case where weighting stops being cosmetic.
One rerun will not show you that. Across draws the unweighted estimate carries the larger bias, but within any single draw the separation is small next to that draw’s own sampling variation, so one comparison is an unreliable demonstration and can as easily show the weighted estimate sitting further away. Run eight or ten seeds, keep each $SIM_OUT file, and fit the anchor to every one of them weighted and unweighted. $SIM_SWEEP will not do this job for you, because it records the unweighted estimates only.
Two things to expect. The weighted estimate will not land exactly on the planted value either: w1 is the reciprocal of the selection probability, which reverses the distortion in expectation rather than in any one sample. And a distortion that varies across cells changes what the per-cell threshold solve was calibrated against, so the cell means in the twin drift slightly from those in the panel – visible in the generator’s own diagnostic, and a reminder that a generator’s parts are not independent of one another.
The inference reading argues that to test a specification curve you must know what the curve would look like if there were nothing to find, and the test designed there builds that nothing by permutation. Simulation offers a more direct route: write the recipe with the effect removed.
# One null world: the cell structure of the twin with the slope written out of
# the recipe, and a small menu of estimators fitted to it. Every estimate below
# is of a quantity that is exactly zero.
simulate_null <- function(seed) {
set.seed(seed)
cells <- panel |>
select(cntry, year, n_cy, unemp) |>
mutate(n_sim = pmax(1L, as.integer(round(n_cy * 0.10))),
u_cy = rnorm(n(), 0, SD_CY))
cells |>
left_join(cells |> distinct(cntry) |>
mutate(u_c = rnorm(n(), 0, SD_CNTRY)), by = "cntry") |>
uncount(n_sim) |>
mutate(cosmo = pmin(pmax(MEAN_COSMO + 0 * unemp + u_c + u_cy +
rnorm(n(), 0, SD_PERSON), 0), 1))
}
d_null <- simulate_null(2026)
null_fits <- list(
pooled = \(d) lm(cosmo ~ unemp, data = d),
country_re = \(d) lmer(cosmo ~ unemp + (1 | cntry), data = d),
three_level = \(d) lmer(cosmo ~ unemp + (1 | cntry) + (1 | cntry:year), data = d),
cell_means = \(d) lm(mcosmo ~ unemp,
data = summarise(d, mcosmo = mean(cosmo),
unemp = first(unemp),
.by = c(cntry, year)))
) |>
imap(\(f, nm) {
p <- parameters::model_parameters(f(d_null))
tibble(estimator = nm, estimate = p$Coefficient[p$Parameter == "unemp"],
p_value = p$p[p$Parameter == "unemp"])
}) |>
list_rbind()
null_fits |> mutate(across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 4 × 3
estimator estimate p_value
<chr> <dbl> <dbl>
1 pooled -0.00132 0.000395
2 country_re 0.00131 0.00914
3 three_level 0.00104 0.188
4 cell_means -0.000988 0.249 Every row estimates a quantity known to be zero, and on this seed 2 of the four reject it at the five per cent level. Rerun with a different seed (the function takes one for exactly that reason) and the two extremes hold. Over thirty draws the pooled row, which treats forty thousand echoing respondents as independent, rejects in about six draws in ten. The three-level row, which counts the units that actually carry information about unemployment, rejects rarely – blocks of thirty draws put it anywhere from none to one in ten – and is the closest of the four to the five per cent that a valid test should deliver. The two rows in between land between those extremes and move about from one block of seeds to the next.
A miniature multiverse fitted to a null world shows what a full-scale permutation test constructs: the distribution of whole curves under nothing. It also disciplines how the percentages of the real curve should be read, a point that the multiverse module develops – some share of any universe will look supportive under a true zero, so a share alone, with no null to stand it against, is a description rather than a verdict.
real_panel <- sample(1:9, 1)
map(1:9, \(i) {
d <- if (i == real_panel) sim_cells else
sim_cells |> mutate(unemp = sample(unemp), .by = cntry)
d |> mutate(panel = i)
}) |>
list_rbind() |>
ggplot(aes(unemp, cosmo)) +
geom_point(alpha = 0.4, size = 0.8, colour = "#00205b") +
geom_smooth(method = "lm", se = FALSE, linewidth = 0.6, colour = "#c8102e") +
facet_wrap(~panel) +
labs(x = "unemployment (%)", y = "cell mean, simulated cosmopolitan share") +
theme_minimal(base_size = 12)`geom_smooth()` using formula = 'y ~ x'
Panel 8. If you picked it, you have informally rejected the null by eye – the logic of a permutation test, run on your visual system. If you could not, that is itself evidence about how strong a level-2 relationship of this size looks at this sample size.
Power analysis for multilevel designs rarely has clean formulas, and the standard advice is to simulate: generate data from the design you are contemplating and see how precisely it estimates the quantity you care about (Gelman et al. 2020; DeBruine and Barr 2021). Precision on a level-2 slope is bought mainly with countries and years rather than with respondents, so the question to ask is what four times as many of each would buy.
simulate_design <- function(k, seed, per_cell = 1) {
set.seed(seed)
cells <- panel |>
select(cntry, year, n_cy, unemp) |>
crossing(clone = seq_len(k)) |>
mutate(cntry = paste0(cntry, "_", clone),
n_sim = pmax(1L, as.integer(round(n_cy * 0.10 * per_cell))),
u_cy = rnorm(n(), 0, SD_CY))
cells |>
left_join(cells |> distinct(cntry) |>
mutate(u_c = rnorm(n(), 0, SD_CNTRY)), by = "cntry") |>
uncount(n_sim) |>
mutate(cosmo = pmin(pmax(MEAN_COSMO + B_COSMO * unemp + u_c + u_cy +
rnorm(n(), 0, SD_PERSON), 0), 1))
}
list(c(k = 1, per_cell = 1), c(k = 1, per_cell = 4), c(k = 4, per_cell = 1)) |>
map(\(d) {
m <- lmer(cosmo ~ unemp + (1 | cntry) + (1 | cntry:year),
data = simulate_design(d[["k"]], seed = 20260727 + d[["k"]],
per_cell = d[["per_cell"]]))
p <- parameters::model_parameters(m, ci_method = "satterthwaite")
tibble(countries = 27 * d[["k"]], respondents = nobs(m),
se = p$SE[p$Parameter == "unemp"],
df = p$df_error[p$Parameter == "unemp"])
}) |>
list_rbind() |>
mutate(across(where(is.numeric), \(x) signif(x, 3)))# A tibble: 3 × 4
countries respondents se df
<dbl> <dbl> <dbl> <dbl>
1 27 41700 0.000895 256
2 27 167000 0.000757 253
3 108 167000 0.000409 1060Read the three rows against each other. Quadrupling the respondents inside the existing country-year map does buy precision (people are not worthless), but it buys much less than quadrupling the countries, and it buys almost no degrees of freedom, because the number of units carrying independent information about unemployment has not changed. No feasible amount of extra survey data inside the existing 27 countries would do what 81 more countries would, and no choice on the specification menu can manufacture level-2 information that the design does not contain. That is why this literature argues within wide uncertainty bands, and simulation measures a design’s ceiling before anyone commits to it. A methodological literature exists on doing this properly – declaring aims, data-generating mechanisms, estimands, methods and performance measures rather than tinkering (Morris et al. 2019) – and it is the natural next reading.
Four things to read next. For the habit itself, Regression and Other Stories is the anchor text, its fake-data simulations short, constant and cumulative (Gelman et al. 2020), with Active Statistics packaging the same moves as classroom activities (Gelman and Vehtari 2024). For multilevel structures specifically, DeBruine and Barr’s tutorial is sections 2 and 6 of this module done properly: it builds a mixed-effects simulation from nothing and recovers its parameters, on the argument that constructing the data is how you come to understand the model (DeBruine and Barr 2021). For simulation inside a modelling workflow – prior and posterior predictive checking, calibration, experimentation on synthetic data before touching the real thing – the Bayesian workflow book is the current statement (Gelman et al. 2026), and its habits transfer well beyond Bayesian settings. And for simulation as a research method in its own right, the design framework of Morris et al. (2019) supplies the vocabulary, the book-length treatment written for social scientists is Carsey and Harden (2014), and an applied Monte Carlo study that this site already cites, on cross-lagged panel models, shows the genre working on a question adjacent to this case’s own estimator fork (Leszczensky and Wolbring 2022).
Simulation studies are analyses too, and they have researcher degrees of freedom of their own – conditions, performance measures, which results get reported – so a comparative study can be tuned until any method looks superior (Pawel et al. 2024). The disciplines that this course keeps returning to, declared designs and justified choices, apply to simulated data exactly as they apply to real data. The generator, the calibration script, the codebook and this page’s chunks are all committed, so the first Try-it is the quickest place to begin: change what the recipe plants, rerun, and watch what your tools do when you already know the answer.