REPLICATION CODE: 1963 DEMONSTRATION
"Beyond Isolated Headlines: A Framework for Coupling Scorekeeping and
Explanatory Modeling in Public Political Polls"

PURPOSE
  Reproduces every table and figure in the paper that uses the 1963
  Harris-Newsweek survey, plus the Appendix B.7 bootstrap and common-basis
  increment comparisons, all three routes of
  the Appendix E.2 triangulation example (model, propensity, and Bayesian
  network), the Auto-Build search that produced the Appendix E.3 alternative
  specification, and the Nonresponse Test, Variable-Ranking Tool, and
  synthetic-variable figures reported in Appendix E.4. The script is standalone: it does not use
  the companion site's software.

FILES
  replicate_1963.py   The replication script.
  opster6.csv         Analysis file: the Harris 1963 Newsweek Presidential
                      Election Survey, Study No. 1285 (UNC Dataverse,
                      hdl:1902.29/H-1285), 1,283 respondents. Columns q1-q42
                      are the archived items (Supplementary Appendix S1). The
                      remaining columns are author-created; see VARIABLES
                      ADDED BY THE AUTHOR below.
  jfk_codebook.csv    Codebook the companion site uses for the 1963 file: one
                      row per variable level. The script reads three columns
                      for the Appendix E.3 search: variable, actionability
                      (the tier, 0-3; the Actionable path keeps tier 3), and
                      type (three weight columns are flagged continuous). The
                      other columns (label, category, bn_tier, level,
                      level_label, prechecked, non_substantive, scale) are the
                      site's display and network metadata and are not used.
  requirements.txt    Python packages.

HOW TO RUN
  pip install -r requirements.txt
  pip install --no-deps pgmpy==0.1.23
  python replicate_1963.py --data opster6.csv --codebook jfk_codebook.csv --out output --check

  Quick verification (a few seconds; see the --core-check paragraph below):

  python replicate_1963.py --data opster6.csv --codebook jfk_codebook.csv --core-check

  pgmpy 0.1.23 is the Bayesian-network library release the companion site
  uses. Its package metadata declares torch, which it never imports, so it is
  installed without dependencies; requirements.txt lists what it does use.
  The script restores the numpy aliases that release expects.

  Run both install lines. Installing pgmpy with --no-deps and skipping
  requirements.txt leaves pgmpy importable but the packages it uses missing.

  If pgmpy or any package it uses is missing, the script skips only the network
  route: it names the missing package, marks the nine network values SKIP, and
  still runs and checks everything else. The notice distinguishes the two cases,
  so "pgmpy is installed but the module opt_einsum it needs is missing" means
  run pip install -r requirements.txt, not reinstall pgmpy.

  --check compares 130 reported values with the computed ones at the paper's
  rounding. Exit status 0 means every value was reproduced; 1 means at least
  one differs; 2 means every value checked was reproduced but the network route
  was unavailable, for the reason the notice gives.

  --core-check is a quick verification, not a replacement for --check.
  Reviewers who want to verify the principal manuscript-facing results
  without running the computationally intensive sensitivity analyses can use
  it: in a few seconds it recomputes and compares 78 core values - the
  reference-model sample and 56.6 percent baseline, the principal Table 1
  scenario results and intervals, the Table 2 standalone-incremental
  comparisons, the Table B1 diagnostics including the cross-validated Tjur
  R-squared, the Table C1 civil-rights coefficient, the Republican subgroup
  result the walkthrough uses, the Table E1 sample and the Appendix E.3
  common-sample comparison, the Appendix A.2 weight and race-raking
  quantities, the Appendix E.2 model- and propensity-route scenario results
  with balance, the nonresponse tipping calculations, and the propensity
  support rule. It does not run or verify the Appendix B.7 bootstraps, the
  Bayesian-network analyses, the synthetic-variable seeds, the full
  variable-ranking analysis, or the Appendix E.3 Auto-Build search; its
  output names those exclusions and defers to --check, which remains the
  authoritative 130-value test. The core values are a tagged subset of the
  script's single expected-values table, so the two checks compare against
  the same published entries and cannot drift apart. --core-check writes no
  output files, never imports the network libraries (so it runs without
  pgmpy), and exits 0 when every core value is reproduced and 1 otherwise.

  Run time is about two and a half minutes, most of it the Appendix B.7 and
  Appendix E.2 bootstraps; the Appendix E.3 search takes about 20 seconds.

OUTPUT (in the --out folder)
  table_1.csv           Table 1
  table_2.csv           Table 2, plus the party-identification and
                        ideology values given in its note
  table_b1.csv          Table B1
  table_b2.csv          Table B2
  table_c1.csv          Table C1
  table_d1.csv          Table D1, plus the Appendix D.2 race margins
  table_e1.csv          Table E1
  sample_summaries.csv  Sample sizes and baselines for Tables 1, D1, and E1,
                        the Appendix E.3 comparison on the 1,147
                        respondents complete on both specifications, and the
                        Appendix A.2 weight summary (range, weighted n, design
                        effect, Kish effective sample size, precision
                        benchmark), and the Appendix B.7 pinned sample size,
                        seed, and reversal share
  b7_common_basis.csv   Appendix B.7 common basis: drop-one increments for
                        q20a and q31_13 from the reference model plus q31_13
                        (1,156 cases), and for q31_9, q21a, and Racenew from
                        the reference model plus q31_9 (1,153 cases), with
                        2.5th and 97.5th bootstrap percentiles for the latter
                        (the shares of resamples in which q31_9 exceeds q21a
                        and Racenew are in sample_summaries.csv)
  b7_bootstrap.csv      Appendix B.7: the four quantities on the pinned
                        1,156-case sample, with 2.5th and 97.5th bootstrap
                        percentiles (the reversal share is in
                        sample_summaries.csv)
  a2_race_raking.csv    Appendix A.2 race-raking sensitivity: model-sample
                        approval, design effect, and civil-rights "Oppose"
                        coefficient under the published weight and under the
                        two race benchmarks
  e2_propensity.csv     Appendix E.2 (and the "Four Vulnerabilities" and
                        walkthrough sections): propensity side of the
                        triangulation example, Civil Rights "Oppose" -
                        adjusted approval, change from the baseline with its
                        bootstrap 95% CI, the paired gap with the model route
                        and its CI, matching respondents and respondents in
                        each propensity quintile, weighted quintile shares,
                        balance after subclassification, and the E-value
  e2_network.csv        Appendix E.2 (and the "Four Vulnerabilities" and
                        walkthrough sections): network side of the
                        triangulation example - scenario and baseline
                        approval, change from the observed baseline, the
                        paired 95% Dirichlet band, and which reference
                        predictors reach approval directly, only through
                        another predictor, or not at all
  e2_network_edges.csv  The learned network's edges with their bootstrap
                        confidence (share of 40 weighted refits recovering
                        the same adjacency in either direction)
  e4_nonresponse.csv    Appendix E.4 Nonresponse Test: precision benchmark,
                        frontier floors, tipping rates at 60 and 3 percent
                        response, the kappa extension, the propensity-
                        correlation bias bound, and the within-survey
                        approval gaps shown for scale
  e4_ranking.csv        Appendix E.4 Variable-Ranking Tool: standalone and
                        incremental weighted Tjur R-squared and their sample
                        sizes for all 123 candidates (matches the export in
                        1963-model-materials to its four printed decimals)
  e4_propensity_support.csv
                        Appendix E.4: for each of the 25 single-level
                        scenarios of the published model, respondents who
                        gave the answer in each propensity quintile and
                        whether the scenario meets the five-per-quintile rule
                        (16 do)
  e4_synthetic.csv      Appendix E.4 synthetic-variable probe: for each of the
                        30 starting seeds, the probe's achieved outcome and
                        largest predictor correlations, the refit model's
                        Tjur R-squared, and the civil-rights "Oppose" change
  e3_autobuild_steps.csv
                        Appendix E.3 Auto-Build search: every candidate scored
                        at every forward step, in scoring order - whether it
                        was a random-forest addition, its cross-validated
                        score, whether the collinearity gate or the separation
                        screen skipped it, and whether it was admitted - then
                        one row per backward-pass decision (weakest
                        predictor, its drop-one loss, removed or kept)
  e3_autobuild_pool.csv
                        Appendix E.3: every column of the file with the
                        screen's verdict (in the pool, or the reason it is
                        not), its actionability tier, its standalone weighted
                        Tjur R-squared, and its random-forest rank
                        (the search's pool sizes, forward order, final score,
                        backward-pass result, whether any random-forest
                        addition was admitted, and whether the selection is
                        the same without the additions are in
                        sample_summaries.csv)
  e2_scenario.csv       Appendix E.2 (and the "Four Vulnerabilities"
                        section): model side of the triangulation example,
                        Civil Rights "Oppose" - respondents giving that
                        answer, predicted approval, and the change from the
                        baseline with its 95% CI (the same values as the
                        Table 1 "Oppose" row)
  figure_1.png          Figure 1
  figure_1_points.csv   The marked points in Figure 1 and the change between
                        each pair

METHOD NOTES
  Model: binomial GLM with the analysis weight as a variance weight, rescaled
  to average 1 in each estimation sample, and Huber-White robust standard
  errors with no finite-sample correction (HC0; statsmodels applies none to
  GLM results under any HCx label). Coefficients match Stata's
  logit ... [pw=Weight_Final], robust. Stata multiplies the same robust
  variance by N/(N-1); for the reference model that makes its standard errors
  larger by at most 0.0004.

  Approval: q7 codes 1-2 (Excellent, Pretty good) = 1; codes 3-5 = 0.
  Each predictor enters as indicators, with its lowest code as reference.

  Scenarios: each row sets one predictor to one level for every respondent,
  keeps the other predictors as observed, and averages predicted
  probabilities on the weight. Changes are measured from the model-sample
  baseline. Delta-method intervals for changes include the covariance
  between scenario and baseline.

  Stability flag (†): fewer than 50 respondents at the level, less than
  3 percent of the weighted sample, or a change interval at least
  30 points wide (Appendix B.4).

  Tjur R-squared: class means on the analysis weight.
  AUC: weighted, with ties counted one-half.
  Calibration: deciles cut on the unweighted distribution of predicted
  probabilities; rates within each decile use the weight.
  Cross-validation: 5 folds, scikit-learn KFold(shuffle=True,
  random_state=42), rows in file order.

  Bootstrap (Appendix B.7): rows complete on approval, the weight, the six
  reference predictors, and q31_13 (1,156). On those rows: standalone Tjur
  R-squared for q31_13 and q20a; the q20a increment (reference minus
  reference without q20a); the q31_13 increment (reference plus q31_13 minus
  reference). 1,000 resamples drawn with replacement using
  numpy RandomState(1963); the weight is renormalized to mean 1 in each
  resample. The reversal share is the percentage of resamples in which q31_13
  is higher standalone and q20a is higher incrementally. With 1,000
  resamples its Monte Carlo standard error is about 0.5 points.

  Common basis (Appendix B.7): Table 2 compares add-one increments for
  omitted items with drop-one increments for included predictors. The two
  comparisons the main text draws are recomputed as drop-one losses from one
  model containing both items, on that model's complete cases: reference plus
  q31_13 (1,156; drop q20a, drop q31_13) and reference plus q31_9 (1,153; drop
  q31_9, q21a, Racenew). The second is bootstrapped with 1,000 resamples from
  numpy RandomState(1963), the weight renormalized in each; reported are the
  percentages of resamples in which the q31_9 increment exceeds q21a's (76)
  and Racenew's (99).

  Race raking (Appendix A.2): within the 1,159-case analytic sample, start
  from Weight_Final and alternately rake the six age-by-sex cells (q38) to
  their shares under Weight_Final across all 1,255 respondents with a valid
  weight, and the Black share (Racenew) to 6 percent or 10.2 percent, until
  the largest weight change is below 1e-12. The analytic sample contains only
  White and Black respondents, so the nonwhite target applies to Racenew.
  Neither raked weight is used anywhere else.

  E.2 scenario: the Table 1 routine with q20a = Oppose for every
  respondent; other predictors as observed; change from the model-sample
  baseline with a delta-method interval that includes the scenario-baseline
  covariance. The respondent count is those in the 1,159 who gave that
  answer (356). The function also accepts several items pinned at once.

  E.2 propensity side: Rosenbaum-Rubin (1984) subclassification. A weighted
  logistic regression of giving the scenario answer (q20a = Oppose) on the
  other five reference predictors gives each respondent's propensity score
  (IRLS with a 1e-8 ridge, the companion site's routine). Respondents are cut
  into five subclasses at the unweighted quintiles of that score, tied scores
  staying below a cut. In each subclass the weighted approval rate among
  respondents who gave the scenario answer is multiplied by the subclass's
  share of the full weighted sample, and the five products are added. The
  change is that sum minus the observed weighted baseline. The estimate is
  refused if any subclass holds fewer than five matching respondents. A
  500-draw bootstrap (seed 12345) resamples respondents and refits both the
  propensity route and the model route on each resample, so the interval on
  their gap reflects the shared sample; resamples in which a subclass holds no
  matching respondent are skipped and counted (none here). Balance is the
  weighted standardized mean difference between matching and non-matching
  respondents, averaged within subclasses by subclass share, reported as the
  largest absolute value among each predictor's indicators. The E-value uses
  the ratio of adjusted approval to the baseline.

  E.2 network side: approval (the binary cut) and the six reference
  predictors, as strings, on the 1,159-case sample. Structure by pgmpy's
  hill-climb search on the BIC score, maximum in-degree 3, with every edge out
  of approval and every edge from a predictor to Racenew blacklisted (race sits
  at tier 0, the ratings at tier 2, approval is a sink). Edge confidence from
  40 refits on weighted bootstrap resamples (NumPy default_rng seed 0).
  Conditional probability tables by weighted BDeu estimation, equivalent
  sample size 5, weights rescaled to average 1. The scenario is read by the
  adjusted-prediction rule: each respondent's probability of approval is the
  column of approval's table indexed by the respondent's parent answers, with
  q20a pinned to Oppose, averaged on the weight. The band comes from 300 draws
  (seed 0) of approval's table from its Dirichlet posterior - concentration =
  fitted probability x (5 + weighted respondents at that parent configuration,
  floored at 20) + the BDeu pseudo-count - reading scenario and baseline from
  the same draw and taking percentiles of their difference. The change is
  reported against the observed weighted baseline, as the site reports it.

  E.4 Nonresponse Test: closed-form arithmetic on the model-sample baseline and
  the Kish-based benchmark; no simulation. The within-survey gaps use the same
  1,159-case sample (q3 Democrat minus Republican, q37a Catholic minus
  Protestant, the range across q36 education categories, and the range across
  the six q38 age-by-sex cells).

  E.4 Variable-Ranking Tool: each candidate enters as categorical indicators.
  Standalone fits use the rows complete on approval, the weight, and the
  candidate. Increments are drop-one for the six reference predictors on the
  1,159-case sample and add-one for the 117 omitted candidates on the rows
  complete on the six predictors and the candidate. q40 (three-category race)
  is collinear with Racenew and has no increment, leaving 116. The screened
  subset is the 26 candidates with fewer than 1,000 cases in their incremental
  sample.

  E.4 synthetic-variable probe: the site's generator, reproduced in the script
  (it matches the site's output value for value). Target outcome correlation
  0.30, predictor-overlap ceiling 0.20, three categories, starting seeds 0-29;
  each starting seed tries s, s+7, s+14, s+21 and keeps the admissible probe
  closest to the target. The probe is added to the reference model as a
  categorical predictor; the refit's in-sample weighted Tjur R-squared and its
  "Oppose" scenario change are reported for each seed. The walkthrough check
  confirms that the largest departure from the published -10.4 is under two
  points (1.7).

  Figure 1 marks the same two pairs as the published figure: logit 0 to 1
  (+23.1 points) and logit 5 to 6 (+0.4 points).

  E.3 Auto-Build search (Auto: Actionable path), reimplemented from the
  companion site's routine on the full 1,283-respondent file with the binary
  approval cut as the outcome and no wall-clock budgets. Candidate pool, in
  file order: a column is excluded if it is not in jfk_codebook.csv, is on
  the outcome's declared exclusion list (q7, q8a-q11b, q30a, q8lean, q7_cat,
  JFKvGoldwater, TaxCuts, and the Approval column itself), is more than 40
  percent missing, is constant, correlates with the outcome above |r| = 0.70
  on its numeric codes, or has a pairwise n with the outcome below half the
  median candidate n. Each remaining candidate is scored alone on its own
  complete cases by weighted Tjur R-squared (the standalone fit of E.4). A
  candidate whose |r| with an already-kept candidate of higher standalone fit
  exceeds 0.85 is dropped as a near-duplicate (on this file that removes only
  Racenew against q40, both outside the tier). The Actionable path then keeps
  the codebook's top actionability tier (3), giving 53 candidates.
    Random-forest ranking: a design over the pool (missing = reference level)
  on the rows with a valid outcome and weight; a stratified 70/30 split
  (scikit-learn train_test_split, random_state 42); RandomForestClassifier
  with 300 trees, min_samples_leaf 5, random_state 42, fit with the analysis
  weight; grouped permutation importance = the mean drop in weighted holdout
  AUC over 5 permutations of each variable's columns, permuted in pool order
  from numpy default_rng(42); ranked by that drop, ties by name.
    Forward pass: each step ranks the remaining pool by actionability tier
  then standalone fit, takes the top 12, and appends the 4 highest-ranked
  remaining random-forest candidates not already listed. Candidates are
  scored in that order by the mean over 5 folds (KFold, shuffle, random_state
  42, over the complete cases in file order) of the out-of-fold Tjur
  R-squared, each fold's model a binomial GLM fit with the analysis weight
  rescaled to mean 1 in the training fold and the fold's class means taken
  unweighted (the site's rule; not the pooled, weighted quantity Table B1
  reports). With at least one predictor selected, a candidate whose
  provisional design gives any GVIF^(1/(2df)) above sqrt(3) (VIF 3) is
  skipped. A candidate is admitted only if its score exceeds the running
  best by 0.005 - so the first candidate to clear the gate holds it unless a
  later one beats it by the full margin - and a weighted fit of the proposed
  set converges with every coefficient within 15 in absolute value. The
  pass stops when no candidate clears the gate, or at 10 predictors.
    Backward pass: with 5 folds, repeatedly remove the predictor with the
  smallest drop-one loss while that loss is below 0.005.
    Result: q31_9, q31_1, q31_20, q31_4, q20a, q31_10, q31_3 admitted in
  that order (step 8's best gain, 0.004, is under the gate); the backward
  pass removes nothing (weakest drop-one loss 0.007, q31_3); 1,176 complete
  cases. No random-forest addition was ever admitted, and the script also
  reruns both passes without the additions and records that the selection is
  unchanged, so the one step whose arithmetic depends on the scikit-learn
  release (the forest) does not decide the specification.

VARIABLES ADDED BY THE AUTHOR
  Used by this script
    Weight_Final      Age-by-sex poststratification weight targeted to the
                      estimated 1963 electorate; the analysis weight used
                      throughout the paper (Appendix A.2). Valid for 1,255.
    Racenew           Two-category race from q40: 1 White, 2 Black; q40
                      "Other" and missing are set to missing (Appendix A.2).

  Not used by this script or the paper's tables
    Approval          q7 recoded 1 = Excellent or Pretty good, 0 = Only fair,
                      Poor, or Not sure. Identical to the coding this script
                      applies to q7.
    JFKvGoldwater     Companion-site vote-intention outcome: Kennedy-Goldwater
                      choice including leaners (q8a, then q8b), 1 Kennedy,
                      0 Goldwater, Not sure missing. It also uses the q8b
                      answer for 9 respondents with no q8a answer.
    q8lean            Companion-site three-category version of the leaned
                      choice: 1 Goldwater, 2 Kennedy, 3 Not sure; missing
                      when q8a is missing.
    TaxCuts           Companion-site tax-cut outcome: q30a recoded 1 Favor, 0 Oppose; Not sure is missing.
    Weight_AgeSexVAP  Age-by-sex weight targeted to the 1963 voting-age
                      population rather than the electorate; retained for
                      comparison with Weight_Final.
    Weight_Race       Superseded race-only weight (White 0.98, Black 1.47),
                      the analysis weight before 2026; Weight_Final replaced
                      it.
    Include           Superseded sample flag (1 for 1,180 respondents; the
                      1,159 complete cases plus 21 others). Retired: the
                      published model's complete-case rule now defines the
                      analytic sample. Not used in the paper.

NOT INCLUDED
  The walkthrough's Inside the Model example (about 18 percent) is arithmetic on
  the Table C1 coefficients and is not written out separately. The companion
  site's software itself is not part of the package; the script reimplements
  each estimator whose results the paper reports.

TESTED WITH
  Python 3.12, numpy 2.4, pandas 3.0, scipy 1.17, statsmodels 0.15,
  scikit-learn 1.8, matplotlib 3.10.
