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