---
title: "Grade 4-to-5 Joint Mathematics–Science Benchmark Transitions in Nine Systems: A Descriptive Common-Severity Analysis Using TIMSS 2023 Longitudinal Data"
subtitle: "Supplementary Material — Additional file 1, anonymous supplementary methods, results, and reproducibility record"
format:
  docx:
    toc: true
    number-sections: false
  pdf:
    pdf-engine: tectonic
    documentclass: scrartcl
    classoption:
      - DIV=11
      - numbers=noendperiod
    papersize: letter
    header-includes:
      - \KOMAoption{captions}{tableheading}
      - \usepackage{needspace}
    keep-tex: false
    toc: true
    number-sections: false
bibliography: references.bib
csl: apa.csl
link-citations: true
lang: en
---

This supplement accompanies the anonymous manuscript *Grade 4-to-5 Joint Mathematics–Science Benchmark Transitions in Nine Systems: A Descriptive Common-Severity Analysis Using TIMSS 2023 Longitudinal Data*. It documents the analytic files, joint estimands, matched four-dimensional plausible-value/jackknife repeated replication (PV/JRR) calculations, six standardization specifications, validation checks, robustness analyses, and a descriptive Grade 8-to-9 extension. Every displayed data-driven table is generated directly from the validated canonical objects that also produce `outputs/tables/` and `outputs/diagnostics/`; complete-precision results remain in the cited machine-readable CSV files. BB denotes below both mathematics and science thresholds, MB mathematics below only, SB science below only, and AA at or above both thresholds. The focal estimand is the BB-to-AA transition to joint attainment, not a claim that students had attained previously.

## S1. Data files, systems, samples, and variable map

### S1.1 Official source files and acquisition

The reproducible acquisition script `analysis/01_acquire_data.R` downloads the official TIMSS 2023 Longitudinal Grade 4-to-5 and Grade 8-to-9 R archives, codebooks, official R programs, user guide, longitudinal technical-report addendum, TIMSS standard-error chapter, six official results workbooks, and two sampling/participation workbooks. Source URLs are declared in `analysis/00_config.R`. The acquisition step writes file size, acquisition status, and MD5 checksums to `data/raw/download_manifest.csv`. Raw archives and extracted public-use files are excluded from the distributable code package; users reacquire them from the official TIMSS site [@fishbein2025technical; @fishbein2026user; @vondavier2025longitudinal].

The reduced Grade 4-to-5 analytic file contains 44,325 students in nine systems: Georgia, Italy, Jordan, the Republic of Korea, Kosovo, Montenegro, North Macedonia, Slovenia, and Sweden. The separate Grade 8-to-9 file contains 16,764 students in Jordan, the Republic of Korea, and Sweden. Grade bands were never pooled.

{{< include generated/supplement_sample_summary.md >}}

Participation at each wave was defined as at least one of the two relevant TIMSS process-participation flags being coded `1`; the loader does not treat arbitrary positive values as participation. Fail-fast gates required every nonmissing source code to be in {1, 2}, where `1` denotes participation, and required the derived wave indicators to be binary. The linked files contained both-wave, baseline-only, and follow-up-only patterns but no students classified as participating in neither wave. Achievement files contained every longitudinal plausible-value cell used by the analyses, including for records whose process flags did not indicate both-wave participation. The primary analysis therefore uses the official PV-completed representation; Sections S7.1–S7.2 examine participation restrictions and departures from it.

### S1.2 Variable map and terminal PV index

The loader requires one achievement record and one linked process record per country–student key. It retains the following fields.

| Role | Grade 4–5 field or pattern | Grade 8–9 field or pattern |
|---|---|---|
| System, school, student identifiers | `IDCNTRY`, `IDSCHOOL`, `IDSTUD` | Same |
| Final student weight | `TOTWGT` | Same |
| JRR zone and member | `JKZONE`, `JKREP` | Same |
| Mathematics baseline PV, draw *j* | `ASMMAT1j` | `BSMMAT1j` |
| Mathematics follow-up PV, draw *j* | `ASMMAT2j` | `BSMMAT2j` |
| Science baseline PV, draw *j* | `ASSSCI1j` | `BSSSCI1j` |
| Science follow-up PV, draw *j* | `ASSSCI2j` | `BSSSCI2j` |
| Supplied mathematics change, draw *j* | `ASMCHG0j` | `BSMCHG0j` |
| Supplied science change, draw *j* | `ASSCHG0j` | `BSSCHG0j` |
| 2023 process participation | `ITPARTPT1_23`, `ITPARTPT2_23` | Same |
| 2024 process participation | `ITPARTPT1_24`, `ITPARTPT2_24` | Same |

Here the terminal digit $j=1,\ldots,5$ is the plausible-value index, as shown in User Guide Chapter 2, printed pp. 51–52 and Exhibit 2.5 [@fishbein2026user]. Continuous changes were calculated directly as follow-up minus baseline within matched draw; supplied change variables were retained for audit. Before any estimate was produced, the loader required finite `TOTWGT`, `JKZONE`, and `JKREP`, strictly positive `TOTWGT`, `JKREP` in {0, 1}, and both members represented within every system-specific jackknife zone. The complete mapping is `outputs/diagnostics/variable_map.csv`.

## S2. Estimands and complete PV/JRR covariance derivation

### S2.1 Joint states and transition estimands

For system $s$, PV draw $j$, wave $t$, and common threshold $c$, let $M_{itsj}$ and $S_{itsj}$ denote mathematics and science PVs. The four-state indicator is

$$
Z_{itsj}=\begin{cases}
BB,&M_{itsj}<c,\ S_{itsj}<c,\\
MB,&M_{itsj}<c,\ S_{itsj}\ge c,\\
SB,&M_{itsj}\ge c,\ S_{itsj}<c,\\
AA,&M_{itsj}\ge c,\ S_{itsj}\ge c.
\end{cases}
$$

Define joint transition validity as

$$
V_{isj}=I\{M_{i0sj},S_{i0sj},M_{i1sj},S_{i1sj}\ \text{are all finite}\}.
$$

For origin $a$ and destination $b$, the draw-specific weighted transition implemented by `transition_vector()` is

$$
\widehat\theta_{ab,sj}=
\frac{\sum_i w_{is}V_{isj}I(Z_{i0sj}=a,Z_{i1sj}=b)}
{\sum_i w_{is}V_{isj}I(Z_{i0sj}=a)}.
$$

The primary estimand uses $a=BB$, $b=AA$, and $c=475$. Each origin row has its own denominator restricted to records with jointly valid baseline and follow-up states; the implementation fails if that denominator is nonpositive and verifies that every row sums to one within $10^{-10}$. Additional estimands are baseline and follow-up state shares; $P(AA_1\mid Z_0\ne AA)$; BB persistence and exit; and mean changes among baseline-BB students in mathematics, science, and $\min(M,S)$. These are population descriptions, not individual diagnoses or causal effects.

### S2.2 Matched four-dimensional plausible-value tuples

The TIMSS 2023 Longitudinal Technical Report states that five achievement PVs were imputed simultaneously from country-specific four-dimensional models covering mathematics 2023, mathematics 2024, science 2023, and science 2024 (Section 6, printed p. 32 and pp. 36–37) [@fishbein2025technical]. Together with the terminal-index evidence in the User Guide, this establishes Grade 4-to-5 draw $j$ as

$$
(\mathrm{ASMMAT1}j,\mathrm{ASMMAT2}j,\mathrm{ASSSCI1}j,\mathrm{ASSSCI2}j),
$$

and Grade 8-to-9 draw $j$ as the analogous B-prefix tuple

$$
(\mathrm{BSMMAT1}j,\mathrm{BSMMAT2}j,\mathrm{BSSSCI1}j,\mathrm{BSSSCI2}j).
$$

Except for the explicit same-wave diagnostic in Section S9, all inferential analyses used one of these five matched tuples. Unequal-index combinations were not treated as primary imputations, and PVs were not averaged into person-level scores before classification [@mislevy1992; @wu2005; @rutkowski2010].

### S2.3 JRR sampling covariance within a matched draw

Input validation required finite design variables, positive final weights, `JKREP` in {0, 1}, and both replicate members in every system-zone. For each finite `JKZONE`, the code constructs two replicates. In the selected zone, the retained `JKREP` member receives factor 2 and the other member factor 0; observations outside that zone retain factor 1. If $\widehat{\boldsymbol\theta}_j$ is a named full-weight vector and $\widehat{\boldsymbol\theta}_{jr}$ its value in replicate $r$, the implemented covariance is

$$
\widehat U_j=\frac{1}{2}\sum_r
(\widehat{\boldsymbol\theta}_{jr}-\widehat{\boldsymbol\theta}_j)
(\widehat{\boldsymbol\theta}_{jr}-\widehat{\boldsymbol\theta}_j)^\top.
$$

The divisor is 2, not the number of replicates. Deviations are from the full-sample estimate. Vector estimation retains off-diagonal covariance [@fishbein2024se].

With $M=5$, pooled point, average sampling covariance, and sample between-PV covariance are

$$
\overline{\boldsymbol\theta}=\frac{1}{M}\sum_{j=1}^{M}\widehat{\boldsymbol\theta}_j,
\qquad
\overline U=\frac{1}{M}\sum_{j=1}^{M}U_j,
$$

$$
B=\frac{1}{M-1}\sum_{j=1}^{M}
(\widehat{\boldsymbol\theta}_j-\overline{\boldsymbol\theta})
(\widehat{\boldsymbol\theta}_j-\overline{\boldsymbol\theta})^\top.
$$

Total covariance is

$$
T=\overline U+\left(1+\frac{1}{M}\right)B.
$$

Standard errors are $\sqrt{\operatorname{diag}(T)}$, and displayed 95% intervals are $\widehat\theta\pm1.96SE$ [@rubin1987].

### S2.4 Primary common-severity standardization

Within each draw and system, the baseline-BB risk set defines scaled deficits

$$
d_M=(475-M_0)/100,\qquad d_S=(475-S_0)/100.
$$

The primary joint-attainment surface is

$$
\operatorname{logit}\{P(Y=1\mid d_M,d_S,s)\}
=\beta_{0s}+\beta_{1s}d_M+\beta_{2s}d_S
+\beta_{3s}d_M^2+\beta_{4s}d_S^2+\beta_{5s}d_Md_S,
$$

where $Y=I(Z_1=AA)$. Positive survey weights are divided by their within-fit mean, a constant scaling that leaves weighted score equations unchanged. For numerical stability, prediction evaluates $\operatorname{logit}^{-1}(\eta)$ after clamping $\eta=X\widehat\beta$ to $[-30,30]$. Diagnostics retain the unclamped minimum and maximum and the pooled, own-system, and common-target counts and weighted shares for which clamping occurs.

For each system, weighted 1st and 99th percentiles are calculated separately for baseline mathematics and science in the BB set. Primary common support is the intersection of the nine marginal intervals. If $E_{isj}=1$ indicates that record $i$ lies in the rectangle, its equal-system reference weight is

$$
q_{isj}=\frac{1}{9}\frac{w_{is}E_{isj}}
{\sum_{h:s_h=s}w_{hs}E_{hsj}}.
$$

Each system contributes mass 1/9, regardless of represented population size. The standardized estimate for target system $s$ is

$$
\widehat\theta^{\mathrm{standardized\ common}}_{sj}=
\sum_iq_{ij}\operatorname{logit}^{-1}(X_{ij}\widehat\beta_{sj}).
$$

Models use each target system’s complete BB risk set; integration uses the common target. The result is a modeled descriptive quantity under a synthetic one-system/one-vote severity distribution, not a world-student population or causal estimand. Its interpretation therefore requires an explicit target and adequate positivity/support [@lane1982; @keiding2014; @westreich2010].

### S2.5 Full 168-parameter stacked standardization covariance

The implemented stack is not limited to the 18 primary raw/standardized quantities. Within each matched draw it contains:

- nine `raw_full` quantities shared across specifications;
- for each of the five specifications covering all nine systems, three quantities per system—`raw_support`, `model_own_support`, and `standardized_common` (5 × 9 × 3 = 135); and
- the same three quantities for eight systems under `leave_korea_out` (8 × 3 = 24).

The total is 9 + 135 + 24 = **168 parameters**. The nine raw and nine primary standardized quantities form an exact 18-parameter subvector and covariance submatrix.

A perturbation in one system can alter common support or target weights used for multiple systems and specifications. For every system/zone/member perturbation, the implementation therefore:

1. changes weights only for the selected system and zone;
2. rebuilds affected 1st–99th percentile bounds and the common rectangle;
3. rebuilds strict 4 × 4 cell boundaries, eligibility, and renormalized target when applicable;
4. rebuilds equal-system or pooled represented-BB target weights;
5. refits the perturbed system’s applicable response surfaces;
6. recomputes every standardized quantity affected by the shared target; and
7. adds one-half the outer product of the full 168-element replicate deviation.

Unperturbed system surfaces remain fixed when their weights do not change, but their standardized quantities are recomputed if the shared target changes. Zones with no BB records have exactly zero risk-set deviation. A replicate is never silently dropped: empty support, an invalid target, rank deficiency, nonconvergence, nonfinite values, or covariance failure stops canonical artifact replacement.

The five 168 × 168 draw-level covariance matrices are pooled as in S2.3. For any contrast vector $a$, $\operatorname{Var}(a^\top\widehat{\boldsymbol\theta})=a^\top Ta$. This supplies covariance-aware primary standardized-minus-raw contrasts and all 44 alternative-minus-primary contrasts. Cross-system pairwise standardized contrasts were retained as an unnumbered exploratory reproducibility diagnostic but excluded from the publication tables because multiplicity adjustment does not remove the risk of rank-oriented interpretation.

## S3. Reproduction of official TIMSS-L statistics

The extraction/weighting audit used 24 official source rows spanning Grade 4-to-5 and Grade 8-to-9 means, growth, and available benchmark percentages. Those rows contain 108 scalar point estimates and 102 numeric standard errors; six additional benchmark standard errors are printed as `~` by TIMSS and were recorded as suppressed rather than numeric. All point and numeric-SE checks passed against published rounded values.

{{< include generated/supplement_validation_summary.md >}}

The source-row comparisons are in `outputs/diagnostics/official_result_validation.csv` and `outputs/diagnostics/official_result_validation_long.csv`. The compact audit is `outputs/diagnostics/audit_summary.txt`; the publication validation-summary file is `outputs/tables/table_s16_validation_summary.csv` and is catalogued in Section S10. Reproducing official marginals validates extraction and weights but not the new joint quantities.

## S4. Full four-state transition matrices

Table S1, `outputs/tables/table_s1_full_transition_matrices.csv`, contains 12 system-by-grade-band matrices × 16 origin–destination cells (192 rows). Table S2, `table_s2_state_distributions.csv`, reports 96 baseline/follow-up state shares. Table S3, `table_s3_continuous_changes.csv`, reports 36 mathematics, science, and weakest-domain changes among baseline-BB students.

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S1.** Complete Grade 4-to-5 transition matrices at 475 among below-both (BB), mathematics-below-only (MB), science-below-only (SB), and at/above-both (AA) states. Cells are matched-PV/JRR percentages; each baseline-state row sums to 100% apart from rounding, columns are follow-up states, and darker shading denotes a larger percentage.

![](../outputs/figures/figure_s1_transition_matrices.png){width=100% fig-alt="Nine four-by-four heatmaps show Grade 4-to-5 transitions from baseline rows to follow-up columns for BB, MB, SB, and AA. Every cell prints its percentage and darker blue denotes a larger value; BB persistence is highest in Jordan and Kosovo, while AA persistence is highest in Korea, Slovenia, and Sweden."}

Each matrix row is conditional on its baseline state and sums to 100% apart from floating-point rounding. Point estimates remain in [0%, 100%]. Because intervals use an untransformed normal approximation, a limit for a sparse nonprimary cell can extend slightly outside that range; limits were not silently truncated. The focal BB-to-AA transition is only one cell: complete matrices distinguish BB persistence, one-domain movement, joint attainment, and movement out of AA.

Across Grade 4-to-5 systems, mean weakest-domain change among baseline-BB students was 18.7–50.0 TIMSS points, while BB-to-AA probabilities were 5.5%–34.1%. Continuous and categorical quantities answer different questions.

## S5. Six standardization specifications, decomposition, and overlap

The primary specification uses all nine systems, equal system mass, the intersection of weighted 1st–99th percentile rectangles, and quadratic-interaction surfaces. Five one-factor sensitivity alternatives were internally frozen before their own estimates were examined, after the primary estimate had been obtained; this internal freeze was not an external preregistration:

1. pooled represented-BB target using original retained student weights;
2. eight-system target leaving the Republic of Korea out, with no Korea estimate;
3. 4 × 4 strict common cells requiring unweighted $n\ge5$ and Kish ESS $\ge3$ in every system;
4. additive linear surfaces; and
5. linear-interaction surfaces without squared terms.

Table S4, `table_s4_standardization_specification_decomposition.csv`, reports all 53 applicable specification-system rows. For each it retains `raw_full`, `raw_support`, `model_own_support`, and `standardized_common`, with the exact identity

$$
\begin{aligned}
&\mathrm{standardized\_common}-\mathrm{raw\_full}\\
={}&(\mathrm{raw\_support}-\mathrm{raw\_full})\\
&+(\mathrm{model\_own\_support}-\mathrm{raw\_support})\\
&+(\mathrm{standardized\_common}-\mathrm{model\_own\_support}).
\end{aligned}
$$

The terms are labeled **support restriction**, **in-support model calibration**, and **modeled target distribution**. They are specification-dependent descriptive/model-based components, not causal effects or a pure composition decomposition.

Table S5, `table_s5_standardization_alternative_contrasts.csv`, reports 44 covariance-aware alternative-minus-primary contrasts. Point differences ranged from −1.0 to +3.7 percentage points, and none of the individual 95% intervals excluded zero. This does not establish equivalence: standard errors reached 7.9 points, and the widest strict-cell interval was −12.3 to +18.6 points. The strict-common-cell alternative retained 10–11 of 16 cells and 75.83%–85.89% of target mass before renormalization, but replicate-wise support rebuilding produced wide intervals. It is an instability sensitivity, not a preferred replacement.

Primary overlap diagnostics assess weighted risk-set retention, Kish ESS, reference mass within each target’s empirical convex hull, and reference mass in occupied cells of a 5 × 5 common-support grid.

{{< include generated/supplement_overlap_summary.md >}}

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S2.** Primary common-support diagnostics across five plausible-value draws. The four panels show retained baseline-BB weight, Kish effective sample size (ESS), equal-system reference mass inside each target system's empirical convex hull, and reference mass in 5 × 5 cells containing at least five target-system cases; boxplots summarize the five draws.

![](../outputs/figures/figure_s2_overlap_diagnostics.png){width=100% fig-alt="Four boxplot panels compare nine systems across five plausible-value draws: retained risk-set weight, Kish effective sample size on a log scale, reference mass inside each system's convex hull, and reference mass in five-by-five cells with at least five target cases. Korea has the lowest effective sample size and markedly the lowest occupied-cell coverage."}

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S3.** Five internally frozen one-factor common-severity standardization alternatives minus the primary estimate. Points are system-specific percentage-point contrasts, bars are covariance-aware 95% confidence intervals, the dashed vertical line marks no difference, and the Republic of Korea is absent only from the leave-Korea-out target.

![](../outputs/figures/figure_s3_standardization_sensitivity.png){width=100% fig-alt="Five contrast panels show alternative-minus-primary percentage points by system, plus a sixth explanatory panel. Every covariance-aware 95% interval crosses the dashed zero line; strict-common-cell intervals are much wider than intervals for the target-weighting and response-surface alternatives."}

Korea is the weakest-overlap case: mean baseline-BB $n$ was 162.8, minimum support ESS 88.35, mean hull coverage 92.7%, and five-case-cell coverage 58.0%. Its raw estimate was 32.5% (95% CI 19.9%–45.1%) and standardized estimate 22.7% (12.9%–32.4%). These quantities require joint qualification and do not support ranking. Table S6, `table_s6_standardization_overlap_by_draw.csv`, contains draw-level primary diagnostics; Table S7, `table_s7_standardization_overlap_summary.csv`, contains the nine-system summary.

The unnumbered reproducibility-only standardized-pairwise-contrast CSV retains 36 covariance-aware pairwise standardized contrasts and Holm-adjusted *p* values for auditability. Its exact basename and SHA-256 digest are recorded in Additional file 3's source manifest. The file is deliberately absent from Additional file 2, the numbered supplement, and the table inventory; it must not be used as a league table.

## S6. Alternative thresholds and crossing buffers

The common-threshold grid was 390, 400, 410, 465, 475, and 485. For each threshold, PV/JRR estimation was repeated for BB-to-AA, BB persistence, BB exit, baseline BB share, mathematics change, science change, and weakest-domain change in both grade bands. Table S8 is the 504-row `outputs/tables/table_s8_threshold_sensitivity.csv`. Main Figure 3 displays Grade 4-to-5 BB-to-AA results and is referenced here rather than embedded again.

At 400, Grade 4-to-5 conditional BB-to-AA estimates were 14.0%–47.0%; at 475 they were 5.5%–34.1%. Korea’s mean baseline-BB count was approximately 26 at 400 and its standard error was 19.5 percentage points. Separately, differences between the published rounded official 400 at/above shares were 0–5 percentage points in mathematics and 1–9 in science. These differences are formed from rounded, whole-cohort, domain-specific marginals; the conditional joint transition uses the baseline-below-both denominator and requires follow-up attainment in both domains.

The 450-to-500 buffer produced Grade 4-to-5 estimates of 0.8%–8.0%. A lower buffer required baseline scores below 380 and follow-up scores at/above 420. Each buffer defines a new risk set and event, not an adjusted primary estimate. Table S9, `table_s9_hysteresis_sensitivity.csv`, contains full results. Korea’s lower-buffer risk set averaged only 15.8 students and is not suitable for stable comparison.

## S7. Missing-data, IPW, and MNAR sensitivity analyses

### S7.1 PV-completed, both-wave, and response-IPW targets

Three analyses have distinct targets:

1. **PV-completed primary:** all sampled longitudinal records represented in public achievement files;
2. **both-wave observed:** students indicated as participating at both waves; and
3. **response-IPW:** baseline participants reweighted for observed follow-up response.

For response-IPW, each system and PV draw uses a weighted quadratic-logit model for follow-up participation with centered baseline mathematics, science, squares, and interaction. If $R_i$ indicates follow-up participation and $\widehat\pi_i$ is fitted probability, respondents receive

$$
w_i^*=w_i\frac{\widehat P_w(R=1)}{\max(.05,\widehat\pi_i)},
$$

and nonrespondents receive zero analysis weight [@robins1994; @seaman2013]. Models, stabilization, probabilities, weights, and the PV-completed/both-wave/IPW estimate stack are rebuilt within every JRR replicate. This joint stack supplies covariance-aware method-minus-primary standard errors and intervals in Table S10, the 252-row `table_s10_observed_ipw_sensitivity.csv`.

The diagnostic file `outputs/diagnostics/response_ipw_diagnostics.csv` reports, by system and PV draw, sample counts, weighted response rate, propensity minimum and weighted 1st/5th/50th/95th/99th percentiles and maximum, the fraction at the .05 floor, stabilized IPW-multiplier minimum and corresponding percentiles and maximum, original and reweighted ESS and their ratio, and maximum absolute standardized mean differences across the five non-intercept response-model terms before and after weighting. These checks assess practical positivity and observed-covariate balance; neither model convergence nor finite weights identifies response on unobserved outcomes [@westreich2010]. All fitted models converged and no observed probability reached .05.

For Grade 4-to-5 BB-to-AA, response-IPW estimates differed from the primary by less than 1 percentage point in eight systems; Sweden differed by −1.29 points. These comparisons cannot address unobserved selection.

### S7.2 Pattern-mixture shifts

Three score-shift scenarios were applied to the PV-completed representation:

- subtract $\delta\in\{0,10,20,40\}$ from both follow-up-domain PVs for baseline-only records;
- subtract $\delta$ from both baseline-domain PVs for follow-up-only records; and
- apply both changes under a combined scenario.

Every state, risk set, estimate, JRR covariance, and PV pool was recomputed. Delta zero exactly reproduced primary estimates and standard errors. Under the combined −40 scenario, Grade 4-to-5 BB-to-AA estimates were 5.3%–32.9%; the largest reduction was 2.2 percentage points in Sweden. No system reached half its delta-zero estimate or crossed its delta-zero lower 95% bound within the planned range.

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S4.** Pattern-mixture sensitivity for baseline-only students. Both follow-up-domain plausible values are shifted by 0, −10, −20, or −40 TIMSS points before all states and estimates are recomputed; negative shifts are imposed missing-not-at-random scenarios, not estimated corrections. Vertical bars are 95% confidence intervals.

![](../outputs/figures/figure_s4_mnar_sensitivity.png){width=100% fig-alt="Nine small-multiple curves share a labeled BB-to-AA percentage axis and plot imposed follow-up shifts from minus 40 to zero points per domain for baseline-only students. Curves decline modestly as the shift becomes more negative, with every displayed confidence interval endpoint retained."}

Table S11 is the 1,008-row `table_s11_mnar_pattern_mixture.csv`; Table S12 is the 36-row `table_s12_mnar_tipping_diagnostics.csv`. These are bounded imposed-shift diagnostics, not identified corrections or causal tipping points [@danielsHogan2008; @littleRubin2019].

## S8. Descriptive Grade 8-to-9 extension

Grade 8-to-9 analyses were run separately in three systems, were never pooled with Grades 4-to-5, and are not an independent replication of the nine-system Grade 4-to-5 estimand.

{{< include generated/supplement_grade89_summary.md >}}

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S5.** Separate descriptive Grade 8-to-9 BB-to-AA extension at the 475-point threshold. Points combine five matched plausible-value draws and bars are JRR 95% confidence intervals. The three-system result is not an independent replication of the Grade 4-to-5 analysis, and Jordan science retains its documented cross-wave organization caution.

![](../outputs/figures/figure_s5_grade89_extension.png){width=90% fig-alt="Forest plot of the separate Grade 8-to-9 BB-to-AA extension for Jordan, the Republic of Korea, and Sweden. Estimates are 5.4, 6.8, and 8.1 percent, respectively; Korea has the widest 95% confidence interval, and all interval endpoints are visible."}

Table S13 is `table_s13_grade89_joint_transitions.csv` (48 matrix cells). Domain-specific below-to-at/above estimates were 12.2%, 30.8%, and 25.4% for mathematics and 13.8%, 22.1%, and 18.6% for science in Jordan, Korea, and Sweden. Table S14, `table_s14_grade89_domain_extension.csv`, includes baseline-below shares, persistence, score changes, standard errors, intervals, and cautions. Jordan science is interpreted separately because TIMSS documents changed assessment/instructional organization [@vondavier2025longitudinal].

## S9. Exploratory same-wave unequal-index classification-discordance diagnostic

For wave $t$, all ordered unequal PV-index pairs were used:

$$
D_{ts}=\frac{1}{20}\sum_{j=1}^{5}\sum_{k\ne j}
P_w(Z_{t,k}=AA\mid Z_{t,j}=BB).
$$

Both domains use draw $j$ to define BB and draw $k$ to define AA. The 20 dependent values are averaged arithmetically. This is the sole exception to matched-tuple pooling. Its standard error and interval contain JRR sampling variance for the 20-pair mean only, conditional on the five released PVs; pair standard deviation, minimum, and maximum are separate descriptive quantities and are not a Rubin between-PV component.

```{=latex}
\Needspace{0.75\textheight}
```

**Figure S6.** Actual matched-PV longitudinal BB-to-AA transitions versus baseline and follow-up exploratory same-wave unequal-index classification-discordance means. Same-wave quantities diagnose classification discordance; they are neither longitudinal estimates nor subtractive corrections. Bars for the diagnostic represent JRR sampling-only 95% intervals conditional on the five released PVs.

![](../outputs/figures/figure_s6_measurement_crossing.png){width=100% fig-alt="Grouped forest plot for nine systems with three labeled symbols per system: actual matched-PV longitudinal transition, baseline same-wave unequal-index classification discordance, and follow-up discordance. Every longitudinal estimate exceeds both diagnostics; diagnostic intervals represent JRR sampling variance only, conditional on five released plausible values."}

In Grades 4-to-5, baseline same-wave discordance was 1.4%–6.8% and follow-up discordance 1.8%–11.1%. Follow-up values were approximately 27%–44% of corresponding longitudinal transitions, which exceeded the diagnostic in every system. Plausible values are posterior population draws rather than repeat forms, so neither pair dispersion nor the same-wave mean is a classical measurement-error estimate. Subtraction would be unjustified [@barnett2005]. Table S15 is `table_s15_measurement_crossing_diagnostic.csv`.

## S10. Computational reproducibility and artifact inventory

### S10.1 Dependency-free execution order

From the project root, the complete analysis pipeline is:

```text
Rscript analysis/01_acquire_data.R
Rscript analysis/02_audit_data.R
Rscript analysis/03_primary_transitions.R
Rscript analysis/04_standardization.R
Rscript analysis/05_sensitivity.R
Rscript analysis/06_outputs.R
```

Equivalently, `Rscript analysis/run_all.R` executes those scripts. Acquisition requires internet access; subsequent steps use local official files. The analysis uses base R only, and `analysis/xlsx_extract.py` uses the Python standard library. The exact executable prerequisites and version-check commands are pinned in `BUILD_RECIPE.txt`. The validated toolchain is R 4.6.1, Python 3.14.3, Quarto 1.10.18, Pandoc 3.10 bundled with Quarto, and Tectonic 0.17.0; R must report `capabilities("cairo")` as `TRUE`. The complete document and replacement-package sequence is:

```text
quarto render manuscript/manuscript.qmd --to docx
quarto render manuscript/manuscript.qmd --to latex
quarto render manuscript/manuscript.qmd --to pdf
quarto render manuscript/supplement.qmd --to docx
quarto render manuscript/supplement.qmd --to pdf
python3 tools/finalize_docx.py
python3 tools/finalize_docx.py --verify-only
python3 tools/build_additional_files.py \
  --replace-generated
python3 tools/build_additional_files.py \
  --verify-only
```

The explicit main LaTeX render writes `manuscript/manuscript.tex`; the subsequent main PDF render uses Tectonic, and the supplement PDF render does not retain an additional TeX deliverable. No HTML format is configured or produced.

For a first build in a clean `submission/` directory, the package-builder command omits `--replace-generated`; that flag is reserved for replacing a complete, internally valid generated set. The publication stage creates all eight displayed table fragments from the same in-memory objects as the machine-readable tables, then installs tables, figures, fragments, and validation diagnostics as one rollback-protected set, with the publication-output validation record committed last. The document finalizer rejects a DOCX, PDF, or required main TeX file older than its QMD, bibliography, CSL, finalizer source, included table fragments, or referenced PNG figures. It requires exactly 8 main and 16 supplement source displays, exact nonempty OMML display counts and ordered semantic markers in both DOCX files, source-ordered equations in retained main TeX, and valid PDF signatures/trailers; it also rejects malformed or visibly mangled equations. It removes only canonical passive custom properties and the empty Pandoc comments scaffold, clears identifying and volatile document properties, rejects real comments, tracked revisions, active content, PII, and local paths, and then applies and verifies double spacing outside editable tables, continuous every-line numbering, page-number fields, and exact source-derived figure descriptions.

The package builder freezes one immutable byte snapshot before semantic validation, records every packaged source byte in `provenance/source_snapshot_manifest.csv`, applies deterministic allowlists and archive budgets, rejects active/review DOCX features and identifying or forbidden content, verifies internal and external SHA-256 manifests, and publishes the three Additional Files as one rollback-protected set with the external manifest committed last. Raw and derived student data are deliberately excluded.

### S10.2 Anonymous software record

The validated analysis used exactly R 4.6.1 with base packages, Python 3.14.3 with the standard library, Quarto 1.10.18, its bundled Pandoc 3.10, and Tectonic 0.17.0. Cairo graphics capability was present and is a fail-fast requirement for publication figures. The package builder reconstructs the seven-line anonymous software record from the current executables and requires exact byte equality with the packaged software-version diagnostic; a runtime-version mismatch therefore stops packaging. The machine-specific session-information record is retained outside the anonymous review package. Operating-system build, architecture, numerical-library paths, local user paths, locale, and time-zone details are not needed for anonymous interpretation.

### S10.3 Validation gates and publication inventory

The pipeline stops before replacing publication artifacts if official checks, covariance checks, convergence/rank checks, support or target checks, primary-identity checks, delta-zero identity checks, value-level table round-trips, generated-fragment equality checks, or figure/table artifact checks fail. The validated run reported 13/13 sensitivity checks, 11/11 standardization-producer checks, completion of 7,660 stacked replicate perturbations and 22,980 surface refits without dropping a replicate, a 168-parameter covariance minimum eigenvalue of −4.275395e−18 against −1e−8, and 17/17 publication-output gates.

`outputs/tables/table_manifest.csv` inventories three main and 16 supplemental publication tables. Table S16 is the validation summary and is the final numbered supplemental table. The unnumbered exploratory pairwise-contrast CSV is retained only in Additional file 3 and is absent from the following publication inventory:

{{< include generated/supplement_table_inventory.md >}}

The four main and six supplemental figure pairs, captions, and relative PDF/300-dpi PNG paths are in `outputs/figures/figure_manifest.csv`. Artifact dimensions, sizes, exact file-set checks, declared row counts, role counts, and pass flags are in `figure_artifact_validation.csv`, `table_artifact_validation.csv`, and `publication_output_validation.csv`.

Raw data and regenerated `.rds` files are excluded from distribution. The anonymous reproducibility archive contains the allowlisted `analysis/` and `manuscript/` sources, generated table fragments, nonidentifying diagnostics, machine-readable result tables, figures, official-source checksums, and the source-snapshot manifest, but not `data/raw/`, `data/derived/`, machine-specific session-information diagnostics, author metadata, or local absolute paths.

## References
