Steps analysed: [0, 10, 20, 40, 100, 150]Robustness across the completeness frontier
How coherent results are when more indicators are bought at the cost of fewer regions
In the original analysis of the paper, we demanded that all indicators included in our machine learning analysis cover all the regions that were included in the subjective well-being data provided by the OECD. This rule reduced the indicators that are available in our analysis from 362, all indicators in our combined database, to 135, those available for all regions. Here, I include more indicators and try to gain even finer view on potential hypotheses on subjective well-being.
Based on our completeness rule, one poorly-covered region therefore removes an indicator for all of them, and dropping the worst-covered regions buys indicators back. The published analysis sits at one point on that trade-off. This document asks whether its conclusions depend on where that point was and to explore if other indicators emerge to be important. This analysis is experimental.
Results are read from cached runs of the full pipeline — predictor selection, tuning, importance, ALE and interactions — computed independently at each step. They will differ from the published figures for two reasons. Random forests are stochastic, and ranger and scikit-learn do not produce identical models even from the same seed, so exact agreement is not achievable. Beyond that, every step but the first uses a different dataset than the one the manuscript analysed. What is being compared here is whether conclusions hold, not whether numbers match.
The trade-off
frontier = load_frontier_summary(WAREHOUSE)
marked = frontier[frontier.n_dropped.isin(steps)]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5))
ax1.plot(frontier.n_dropped, frontier.n_indicators, marker="o",
color="black", markersize=4)
ax1.scatter(marked.n_dropped, marked.n_indicators, s=110, facecolor="none",
edgecolor="firebrick", linewidth=2, zorder=3)
ax1.set_ylabel("Indicators available")
ax2.plot(frontier.n_dropped, frontier.n_cells, marker="o",
color="black", markersize=4)
ax2.scatter(marked.n_dropped, marked.n_cells, s=110, facecolor="none",
edgecolor="firebrick", linewidth=2, zorder=3)
ax2.set_ylabel("Total observations (regions x indicators)")
for ax in (ax1, ax2):
ax.set_xlabel("Regions dropped")
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
Dropping the first 10 regions (2.6% of the sample) gains 44 indicators.
Every further step combined gains 45.
Total observations peak at n_dropped = 10.| n_dropped | n_regions | n_indicators | n_cells | |
|---|---|---|---|---|
| 0 | 0 | 388 | 135 | 52380 |
| 1 | 10 | 378 | 179 | 67662 |
| 2 | 20 | 368 | 183 | 67344 |
| 3 | 30 | 358 | 183 | 65514 |
| 4 | 40 | 348 | 183 | 63684 |
| 5 | 50 | 338 | 183 | 61854 |
| 6 | 60 | 328 | 183 | 60024 |
| 7 | 70 | 318 | 183 | 58194 |
| 8 | 80 | 308 | 184 | 56672 |
| 9 | 90 | 298 | 184 | 54832 |
| 10 | 100 | 288 | 199 | 57312 |
| 11 | 110 | 278 | 201 | 55878 |
| 12 | 120 | 268 | 211 | 56548 |
| 13 | 130 | 258 | 216 | 55728 |
| 14 | 140 | 248 | 224 | 55552 |
| 15 | 150 | 238 | 224 | 53312 |
Nearly all of the available gain arrives immediately, and the total quantity of data peaks early. Beyond that point the frontier trades observations away without recovering breadth, so the later steps are useful as a stress test rather than as candidate datasets.
What is traded away
Regions are dropped worst-covered first, and coverage is not missing at random.
worst = dropped_regions(WAREHOUSE, max(steps))
by_country = (
worst.assign(country=worst.reg_id.str[:2])
.country.value_counts()
.rename("regions dropped")
.to_frame()
)
by_country.head(12)| regions dropped | |
|---|---|
| country | |
| ME | 32 |
| TR | 26 |
| US | 24 |
| NZ | 14 |
| CL | 13 |
| JP | 10 |
| AU | 8 |
| KR | 7 |
| IL | 6 |
| EE | 5 |
| IS | 2 |
| FR | 2 |
Where dropped regions cluster in particular countries, the frontier trades indicator breadth for geographic representativeness — a cost the indicator count alone does not show.
Which predictors survive
selected = {n: set(r["selected"]) for n, r in results.items()}
every = sorted(set().union(*selected.values()))
matrix = pd.DataFrame(
{n: [ind in selected[n] for ind in every] for n in steps}, index=every
)
matrix = matrix.loc[
matrix.assign(total=matrix.sum(axis=1))
.sort_values(["total"] + steps, ascending=False)
.index
]
fig, ax = plt.subplots(figsize=(11, 0.26 * len(matrix) + 2))
ax.imshow(matrix.to_numpy(), aspect="auto", cmap="Greys", vmin=0, vmax=1.4)
ax.set_xticks(range(len(steps)))
ax.set_xticklabels(steps)
ax.set_yticks(range(len(matrix)))
ax.set_yticklabels([label(i) for i in matrix.index], fontsize=7)
ax.set_xlabel("Regions dropped")
for edge in range(len(steps) + 1):
ax.axvline(edge - 0.5, color="white", linewidth=2)
plt.tight_layout()
36 distinct indicators selected somewhere; 3 at every step, 9 at exactly one.
Selected at every step:
Disposable income
Sex ratio 65+ (male/female)
Subjective perceived corruptionAn indicator selected everywhere is robust to the trade-off. One appearing once is specific to the indicators that happened to be available there.
The same comparison by concept
Indicator-level churn overstates instability. The richer datasets contain finer-grained variants of the same underlying quantity — EMP_RA gives way to measures such as EMP_SH_PT_F — so an indicator disappearing from the matrix above may mean a more specific member of its family was chosen instead, not that the concept stopped mattering.
concept_selected = {n: {concept(i) for i in r["selected"]} for n, r in results.items()}
every_concept = sorted(set().union(*concept_selected.values()))
cmatrix = pd.DataFrame(
{n: [c in concept_selected[n] for c in every_concept] for n in steps},
index=every_concept,
)
cmatrix = cmatrix.loc[
cmatrix.assign(total=cmatrix.sum(axis=1))
.sort_values(["total"] + steps, ascending=False)
.index
]
fig, ax = plt.subplots(figsize=(11, 0.34 * len(cmatrix) + 2))
ax.imshow(cmatrix.to_numpy(), aspect="auto", cmap="Greys", vmin=0, vmax=1.4)
ax.set_xticks(range(len(steps)))
ax.set_xticklabels(steps)
ax.set_yticks(range(len(cmatrix)))
ax.set_yticklabels(cmatrix.index, fontsize=9)
ax.set_xlabel("Regions dropped")
for edge in range(len(steps) + 1):
ax.axvline(edge - 0.5, color="white", linewidth=2)
plt.tight_layout()
members = pd.DataFrame(
[{"n_dropped": n, "concept": concept(i), "indicator": i}
for n, r in results.items() for i in r["selected"]]
)
(members.groupby(["concept", "n_dropped"])
.indicator.apply(lambda s: ", ".join(sorted(s)))
.unstack()
.fillna("-"))| n_dropped | 0 | 10 | 20 | 40 | 100 | 150 |
|---|---|---|---|---|---|---|
| concept | ||||||
| Age structure | KID_WOM_RA_T, YOU_DEP_RA_T | KID_WOM_RA_T, YOU_DEP_RA_T | KID_WOM_RA_T, YOU_DEP_RA_T | YOU_DEP_RA_T | KID_WOM_RA_T | - |
| Air quality | AIR_POL | AIR_POL | AIR_POL | AIR_POL | - | - |
| Broadband access | BB_ACC | - | - | BB_ACC | BB_ACC | BB_ACC |
| Civic participation | VOTERS_SH | - | - | - | - | - |
| Crime | HOMIC_RA | - | - | - | - | - |
| Education | EDU38_SH | EDU_LF_ISCED58_SH | - | EDU_LF_ISCED58_SH | - | EDU38_SH, EDU_LF_ISCED02_SH |
| Employment | EMP_RA | EMP_RA, EMP_RA_F, EMP_RA_M, EMP_RA_T | EMP_RA, EMP_RA_F, EMP_RA_M, EMP_RA_T | EMP_RA, EMP_RA_F, EMP_RA_T | EMP_PT_SH_F, EMP_PT_SH_M, EMP_PT_SH_T | EMP_PT_SH_F, EMP_PT_SH_M, EMP_PT_SH_SEXDIF_T, ... |
| Housing | ROOMS_PC | ROOMS_PC | - | - | - | ROOMS_PC |
| Income | INCOME_DISP | INCOME_DISP | INCOME_DISP | INCOME_DISP | INCOME_DISP | INCOME_DISP |
| LT_UNEM_RA_T | - | - | - | - | - | LT_UNEM_RA_T |
| LT_UNEM_SH_T | - | - | - | - | - | LT_UNEM_SH_T |
| Mortality | DEATH_RA_M, DEATH_RA_T | DEATH_RA_M | - | - | DEATH_RA_M | - |
| PARTIC_RA_F | - | PARTIC_RA_F | PARTIC_RA_F | PARTIC_RA_F | - | PARTIC_RA_F |
| PARTIC_RA_SEXDIF_T | - | PARTIC_RA_SEXDIF_T | PARTIC_RA_SEXDIF_T | PARTIC_RA_SEXDIF_T | PARTIC_RA_SEXDIF_T | - |
| PARTIC_RA_T | - | PARTIC_RA_T | PARTIC_RA_T | PARTIC_RA_T | - | - |
| Population structure | POP_80M_SH_T, POP_DEN_GR_T, POP_TOT_GI_T | POP_DEN_GR_T, POP_TOT_GI_T | POP_DEN_GR_T, POP_TOT_GI_T | POP_DEN_GR_T, POP_TOT_GI_T | POP_DEN_GR_T | - |
| SURF_T | SURF_T | - | - | - | - | - |
| Sex ratio | SR_ELD_RA_T, SR_TOT_RA_T | SR_ELD_RA_T, SR_TOT_RA_T | SR_ELD_RA_T | SR_ELD_RA_T | SR_ELD_RA_T | SR_ELD_RA_T |
| Subjective measures | SUBJ_PERC_CORR, SUBJ_SOC_SUPP | SUBJ_PERC_CORR, SUBJ_SOC_SUPP | SUBJ_PERC_CORR, SUBJ_SOC_SUPP | SUBJ_PERC_CORR, SUBJ_SOC_SUPP | SUBJ_PERC_CORR | SUBJ_PERC_CORR, SUBJ_SOC_SUPP |
| Unemployment | UNEM_RA | UNEM_RA, UNEM_RA_F | UNEM_RA_F | UNEM_RA_F | UNEM_RA_F | UNEM_RA_F |
Which member of a family is chosen is itself informative: a shift from a general rate to a specific sub-measure says something about which aspect of the concept carries the signal.
Indicators with no concept family assigned (add them to CONCEPT in analysis/labels.py):
LT_UNEM_RA_T
LT_UNEM_SH_T
PARTIC_RA_F
PARTIC_RA_SEXDIF_T
PARTIC_RA_T
SURF_TImportance of the manuscript’s predictors
The manuscript argues that the sex ratio among the elderly is a predictor of comparable importance to disposable income. Tracking that across the frontier tests the claim directly. Importance is aggregated to the concept — the strongest selected member of each family — so a line does not break when the specific measure changes.
imp = pd.concat(
[r["importance"].assign(n_dropped=n) for n, r in results.items()],
ignore_index=True,
)
imp["concept"] = imp.feature.map(concept)
by_concept = imp.groupby(["concept", "n_dropped"]).importance.max().reset_index()
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4.5))
for name in HEADLINE_CONCEPTS:
d = by_concept[by_concept.concept == name].sort_values("n_dropped")
if not d.empty:
ax1.plot(d.n_dropped, d.importance, marker="o", label=name)
ax1.axhline(1.0, color="grey", linestyle=":", linewidth=1)
ax1.set_ylabel("Permutation importance (MSE ratio)")
ax1.set_title("By concept", fontsize=10)
ax1.legend(fontsize=8, frameon=False)
for feature in MANUSCRIPT:
d = imp[imp.feature == feature].sort_values("n_dropped")
if not d.empty:
ax2.plot(d.n_dropped, d.importance, marker="o", label=label(feature, 24))
ax2.axhline(1.0, color="grey", linestyle=":", linewidth=1)
ax2.set_title("The manuscript's specific indicators", fontsize=10)
ax2.legend(fontsize=7, frameon=False)
for ax in (ax1, ax2):
ax.set_xlabel("Regions dropped")
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
by_concept[by_concept.concept.isin(HEADLINE_CONCEPTS)].pivot_table(
index="concept", columns="n_dropped", values="importance"
)| n_dropped | 0 | 10 | 20 | 40 | 100 | 150 |
|---|---|---|---|---|---|---|
| concept | ||||||
| Employment | 2.532 | 1.774 | 1.320 | 1.438 | 1.444 | 1.302 |
| Income | 1.355 | 1.360 | 1.347 | 1.336 | 2.744 | 4.566 |
| Sex ratio | 1.856 | 2.809 | 1.504 | 1.488 | 1.327 | 1.273 |
| Subjective measures | 1.638 | 1.310 | 1.331 | 1.474 | 2.591 | 2.138 |
A gap in the right-hand panel means that specific indicator was not selected at that step. A gap on the left means no member of the family was.
Interaction structure
The manuscript’s central figure rests on the employment rate having the strongest interaction, and on its interaction with the elderly sex ratio. Both the leading predictor and its top partner are re-derived at each step.
rows = []
for n, r in results.items():
top = r["interactions"].head(2)
h1, h2 = float(top.iloc[0].interaction), float(top.iloc[1].interaction)
rows.append({
"n_dropped": n,
"strongest": label(r["top_feature"]),
"H": h1,
"second": label(top.iloc[1].feature),
"H (2nd)": h2,
"margin": h1 - h2,
"margin %": (h1 - h2) / h2 if h2 else np.nan,
"top partner": label(r["top_partner"]),
})
interaction_table = pd.DataFrame(rows).set_index("n_dropped")
interaction_table| strongest | H | second | H (2nd) | margin | margin % | top partner | |
|---|---|---|---|---|---|---|---|
| n_dropped | |||||||
| 0 | Employment rate | 0.242 | Subjective social support | 0.155 | 0.087 | 0.561 | Subjective social support |
| 10 | Sex ratio 65+ (male/female) | 0.227 | Employment rate | 0.119 | 0.109 | 0.916 | Employment rate |
| 20 | Sex ratio 65+ (male/female) | 0.109 | Subjective social support | 0.099 | 0.010 | 0.102 | Employment rate |
| 40 | EDU_LF_ISCED58_SH | 0.100 | Employment rate | 0.091 | 0.009 | 0.097 | Employment rate |
| 100 | Disposable income | 0.160 | Subjective perceived corruption | 0.129 | 0.031 | 0.242 | EMP_PT_SH_T |
| 150 | Subjective social support | 0.178 | Disposable income | 0.157 | 0.021 | 0.137 | Disposable income |
The margin matters more than the ranking. Where two predictors are within a few per cent of one another, which leads is not a stable result — it will move with the seed and with the resolution of the interaction grid.
inter = pd.concat(
[r["interactions"].assign(n_dropped=n) for n, r in results.items()],
ignore_index=True,
)
inter["concept"] = inter.feature.map(concept)
inter_by_concept = inter.groupby(["concept", "n_dropped"]).interaction.max().reset_index()
fig, ax = plt.subplots(figsize=(9, 5))
for name in HEADLINE_CONCEPTS:
d = inter_by_concept[inter_by_concept.concept == name].sort_values("n_dropped")
if not d.empty:
ax.plot(d.n_dropped, d.interaction, marker="o", label=name)
ax.set_xlabel("Regions dropped")
ax.set_ylabel("Interaction strength (H)")
ax.legend(fontsize=8, frameon=False)
ax.spines[["top", "right"]].set_visible(False)
plt.tight_layout()
The interaction surface
The manuscript’s key figure is the joint ALE of the two strongest interacting predictors. Because both are re-derived at each step, the pair being plotted can differ — which is itself the test of whether the published pairing was an artefact of the indicators available.
fig, axes = plt.subplots(1, len(steps), figsize=(4.2 * len(steps), 4))
axes = np.atleast_1d(axes)
for ax, n in zip(axes, steps):
r = results[n]
s = r["ale_2d"]
mesh = ax.pcolormesh(s["x1"], s["x2"], s["values"].T, cmap="Greys", shading="auto")
fig.colorbar(mesh, ax=ax, fraction=0.046)
ax.set_title(f"dropped {n}", fontsize=9)
ax.set_xlabel(label(r["top_feature"], 22), fontsize=7)
ax.set_ylabel(label(r["top_partner"], 22), fontsize=7)
ax.tick_params(labelsize=6)
plt.tight_layout()
dropped 0: Employment rate x Subjective social support - 173/400 grid cells interpolated
dropped 10: Sex ratio 65+ (male/female) x Employment rate - 164/400 grid cells interpolated
dropped 20: Sex ratio 65+ (male/female) x Employment rate - 164/400 grid cells interpolated
dropped 40: EDU_LF_ISCED58_SH x Employment rate - 181/400 grid cells interpolated
dropped 100: Disposable income x EMP_PT_SH_T - 232/400 grid cells interpolated
dropped 150: Subjective social support x Disposable income - 227/400 grid cells interpolatedA high share of interpolated cells means the surface is largely smoothed from neighbouring values rather than estimated from observations. Structure appearing only in heavily interpolated regions should not be read as a finding.
Functional form
fig, axes = plt.subplots(1, len(MANUSCRIPT), figsize=(4.2 * len(MANUSCRIPT), 3.8))
axes = np.atleast_1d(axes)
for ax, feature in zip(axes, MANUSCRIPT):
for n, r in results.items():
d = r["ale_1d"]
d = d[d.feature == feature]
if not d.empty:
ax.step(d.x, d.ale, where="post", linewidth=1.3, alpha=0.85, label=str(n))
ax.axhline(0, color="grey", linestyle=":", linewidth=0.8)
ax.set_title(label(feature, 26), fontsize=8)
ax.tick_params(labelsize=7)
ax.spines[["top", "right"]].set_visible(False)
axes[0].set_ylabel("ALE of subjective well-being")
axes[-1].legend(title="dropped", fontsize=7, frameon=False)
plt.tight_layout()
Curves keeping their shape across steps describe a relationship the model finds regardless of which other indicators are available. Curves changing sign or shape do not.
What this shows
Concepts selected at every step (5):
Employment
Income
Sex ratio
Subjective measures
Unemployment
Of the manuscript's 4 predictor families, 4 survive every step.
At 1 of 6 steps the strongest and second strongest interactions are within 10% of one another.Three things follow.
The published dataset sat one step from a substantially richer one. Dropping a small fraction of the worst-covered regions buys a large share of the available indicators, and total data peaks there rather than at zero. That is a design choice the original analysis made implicitly, by applying a completeness rule to no introduce a selection bias on the covered regions.
Concept-level stability is the meaningful test, not indicator-level. The richer datasets substitute finer-grained measures within the same families, so raw selection churn overstates instability. Read the concept matrix and the members table together.
Interaction rankings are less separated than a single run suggests. Where the margin between first and second place is small, the identity of the strongest interaction is not established, and conclusions resting on that ranking should be stated with the margin attached.
NoteHow to read these comparisons
The datasets are nested subsets of one another, the same data drives both selection and evaluation, and the retained sample changes with every step. These are exploratory robustness checks, not hypothesis tests. A predictor selected at every step is more credible than one selected once; a difference between two steps is not a measured effect.