Use this Skill to analyze PISA/TIMSS international large-scale assessment data: plausible values averaging, Fay's BRR variance estimation, and cross-national SES gradients.
Install with Codex or Claude Copy this prompt, paste it into Codex, Claude, or another assistant, and let it review the skill page and install it for you.
A direct command skips the review prompt. Inspect the source before running it.
Use this Skill to analyze PISA/TIMSS international large-scale assessment data: plausible values averaging, Fay's BRR variance estimation, and cross-national SES gradients.
Do NOT use naive averaging of plausible values without BRR weights —
this underestimates standard errors and produces invalid significance tests.
Background
Why Plausible Values Exist
PISA does not administer every item to every student (matrix sampling design).
Each student answers only a booklet subset, so the full proficiency scale
cannot be directly estimated. Instead, PISA generates M = 10 plausible
values (PV1MATH–PV10MATH) per student — random draws from the posterior
distribution of proficiency given the student's responses and background.
Correct Estimation Protocol (Rubin's Rules)
For statistic Q (e.g., country mean):
1. For each plausible value m = 1..M:
Q_m = estimate using PV_m with BRR replicate weights → Var_BRR_m
2. Point estimate: Q = (1/M) * Σ Q_m
3. Sampling variance: U = (1/M) * Σ Var_BRR_m
4. Imputation variance: B = (1/(M-1)) * Σ (Q_m - Q)²
5. Total variance: Var_total = U + (1 + 1/M) * B
6. SE = sqrt(Var_total)
Fay's BRR Variance Estimation
PISA provides 80 balanced repeated replication (BRR) weights (W_FSTR1–W_FSTR80).
For each replicate r:
defexample_country_means():
"""Compute PISA math means with correct BRR SE for 10 countries."""
df = load_pisa_sample(n=5000)
means = pisa_mean_se(df, group_col="CNT", weight_col="SENWT")
print("PISA Math Means (PV-correct BRR SE):")
print(means.to_string(index=False))
plot_cross_national_comparison(means, output_path="pisa_math_means.png")
return means
if __name__ == "__main__":
example_country_means()
Example 2 — SES Gradient per Country
defexample_ses_gradients():
"""Estimate SES gradient slopes across countries and visualise."""
df = load_pisa_sample(n=8000)
gradients = ses_gradient_by_country(df)
print("SES Gradient Slopes (score points per ESCS unit):")
print(gradients.to_string(index=False))
# Plot slope comparison
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(gradients["group"], gradients["slope"],
xerr=gradients["se"] * 1.96, capsize=4,
color=["#d62728"if s > 40else"#2ca02c"for s in gradients["slope"]])
ax.axvline(0, color="black", linewidth=0.8)
ax.set_xlabel("SES Gradient Slope (score points / ESCS)")
ax.set_title("PISA Math SES Gradient by Country")
ax.grid(axis="x", alpha=0.3)
plt.tight_layout()
plt.savefig("ses_gradients.png", dpi=150)
return gradients
if __name__ == "__main__":
example_ses_gradients()
Example 3 — Gender Gap Analysis
defexample_gender_gap():
"""Compute gender gap in math for each country using PV-correct estimates."""
df = load_pisa_sample(n=6000)
# Code gender: 1=female, 2=male (PISA convention)
results = []
for cnt in df["CNT"].unique():
sub = df[df["CNT"] == cnt]
for gender, label in [(1, "Female"), (2, "Male")]:
sg = sub[sub["ST004D01T"] == gender]
pv_means = [np.average(sg[pv], weights=sg["SENWT"]) for pv in PV_COLS]
results.append({"country": cnt, "gender": label, "mean": np.mean(pv_means)})
gap_df = pd.DataFrame(results).pivot(index="country", columns="gender", values="mean")
gap_df["gap_M_minus_F"] = gap_df["Male"] - gap_df["Female"]
gap_df = gap_df.sort_values("gap_M_minus_F", ascending=False)
print("Gender gap in PISA math (M - F, score points):")
print(gap_df["gap_M_minus_F"].to_string())
return gap_df
if __name__ == "__main__":
example_gender_gap()
