WP35 Phase 3: formal framework (5 theorems, 3 conjectures), D/K/P indices with IRT alternative, measurement diagnostics, a-priori power

Co-Authored-By: Claude Fable 5.1 <noreply@anthropic.com>

Simon-Pierre Boucher committed yesterday (Nov 1, 2026) parent 8621449 Browse files

16 changed files +1,139 −0

added data/processed/modern_indices.parquet +0 −0

Binary file not shown.

added data/processed/seshat_crisis_indices.parquet +0 −0

Binary file not shown.

added data/processed/seshat_indices.parquet +0 −0

Binary file not shown.

added paper/figures/fig_irt_invariance.pdf +0 −0

Binary file not shown.

added paper/figures/fig_measurement_DK.pdf +0 −0

Binary file not shown.

modified research/bib_overrides.csv +2 −0

		`@@ -47,3 +47,5 @@ marzen2017\|10.1098/rsif.2017.0166\|Marzen & DeDeo 2017 J R Soc Interface — adde`
47	47	flack2012\|10.1098/rstb.2011.0214\|Flack 2012 Phil Trans B — added after novelty search
48	48	hoyer2024\|10.31235/osf.io/rk4gd\|Hoyer et al. 2024 SocArXiv preprint 'All crises are unhappy in their own way' — no journal version found 2026-11-01
49	49	wolpert2017\|10.1017/9781316584200.010\|Wolpert 2017 'The many faces of state space compression' (same volume)
	50	+lindley1952\|10.1017/s0305004100027638\|Lindley 1952 Math. Proc. Cambridge Phil. Soc. 48 — verified by Crossref (used in framework.md Theorem 5)
	51	+loynes1962\|10.1017/s0305004100036781\|Loynes 1962 Math. Proc. Cambridge Phil. Soc. 58 — verified by Crossref

added research/framework.md +264 −0

		`@@ -0,0 +1,264 @@`
	1	+# WP35 — Formal framework: coordination demand, coordination capacity, and institutions as codes
	2	+
	3	+Status labels used throughout: DEFINITION (a convention), THEOREM / PROPOSITION (a statement proved below
	4	+or a standard result with the proof reproduced), CONJECTURE (a statement we believe is worth testing, with no
	5	+proof; what would falsify it is stated). Nothing in this file is an empirical finding. The results concern an
	6	+idealized communication system; §8 states what is not claimed and which pre-registered test (if any) can speak to
	7	+each statement. Citations use the keys of `research/reference_candidates.csv`.
	8	+
	9	+---
	10	+
	11	+## 1. Primitives and definitions
	12	+
	13	+DEFINITION 1.1 (agents, tasks, plans). A society has $N$ agents. In each period, coordination tasks arrive at rate
	14	+$\lambda$ per agent. A task requires that $k \ge 2$ participants settle on a common element of a finite plan set
	15	+$S$, $\|S\| = s$, the required plan being drawn independently across tasks from a distribution $p$ on $S$ with Shannon
	16	+entropy $H(p) = -\sum_{x \in S} p(x)\log_2 p(x) \le \log_2 s$ (shannon1948; cover2006). Communication is counted in
	17	+messages (one transmission between two agents) and in bits (the length of what is transmitted). Verification and
	18	+record-keeping use a second channel whose cost is counted separately as administrative overhead (DEFINITION 2.1).
	19	+
	20	+DEFINITION 1.2 (coordination demand). The coordination demand of a period is the task volume that must be
	21	+settled: $D = N \lambda k$ participant-slots per period. In the worst case in which every pair of agents must be
	22	+reconciled, $D$ is bounded by the number of pairs, $N(N-1)/2$, i.e. $D = \Theta(N^2)$; under bounded $\lambda$ and
	23	+$k$, $D = \Theta(N)$. Either way $D$ is increasing in $N$.
	24	+
	25	+DEFINITION 1.3 (coordination capacity). Fix a tolerated error rate $\varepsilon \in (0,1)$. The coordination
	26	+capacity $K$ of an institutional system is the maximal number of tasks it can settle per period such that the
	27	+fraction of tasks settled on a wrong or inconsistent plan is at most $\varepsilon$. $K$ depends on the available
	28	+communication budget (bits per period), the code in use (§2), the routing structure (§5) and the overhead (§3).
	29	+
	30	+DEFINITION 1.4 (utilization). $\rho = D/K$. When $D$ and $K$ are measured in the same units (tasks per period),
	31	+$\log D - \log K = \log \rho$ exactly.
	32	+
	33	+Limitation (not a theorem). The pre-registered empirical index $P = z(\log D) - z(K)$ uses standardized proxies
	34	+(population, territory, settlement hierarchy for $D$; shares of administrative, legal/monetary, communication and
	35	+record-keeping indicators for $K$). It is a monotone proxy for $\log\rho$ only under CONJECTURE C1 (§7); it is not a
	36	+measurement of $\rho$, and the Seshat/V-Dem indicators are not Shannon quantities.
	37	+
	38	+---
	39	+
	40	+## 2. Theorem 1 — coding gain of a shared code
	41	+
	42	+DEFINITION 2.1 (code). A shared code is an injective map $c: S \to \{0,1\}^*$ known to all participants,
	43	+with lengths $\ell(x) = \|c(x)\|$. It is prefix-free if no codeword is a prefix of another. Institutions in the sense of
	44	+this paper (standardized names, units, legal forms, accounting categories, calendars) are shared codes in exactly this
	45	+sense (arrow1974; cremer2007).
	46	+
	47	+Assumption 2.2. Settling a task means that the initiating participant identifies the required plan $x \in S$ to the
	48	+other $k-1$ participants without ambiguity; messages to different participants are counted separately in messages
	49	+but identically in bits per recipient, so we state the result per recipient.
	50	+
	51	+THEOREM 1 (coding gain).
	52	+(i) Without a shared code, identifying an arbitrary element of $S$ without error requires at least
	53	+$\lceil \log_2 s \rceil$ bits per recipient in the worst case.
	54	+(ii) With a shared prefix-free code that is optimal for $p$, the expected length per recipient
	55	+$L_p^\* = \min_{c\ \mathrm{prefix\text{-}free}} \sum_x p(x)\ell(x)$ satisfies
	56	+$$H(p) \le L_p^\* < H(p) + 1 .$$
	57	+(iii) Hence the expected saving per recipient, $G(p) = \lceil\log_2 s\rceil - L_p^\*$, obeys
	58	+$G(p) > \lceil\log_2 s\rceil - H(p) - 1 \ge 0$ whenever $\lceil\log_2 s\rceil \ge H(p)+1$, and $G(p)$ is strictly
	59	+positive as soon as $p$ is sufficiently concentrated ($H(p) < \lceil\log_2 s\rceil - 1$).
	60	+
	61	+Proof. (i) Without a shared code there is no prior agreement on which binary string denotes which plan, so the
	62	+only error-free strategy is a fixed-length enumeration of $S$ agreed by the trivial convention of listing; any scheme
	63	+that distinguishes $s$ alternatives with binary strings of a fixed length $m$ needs $2^m \ge s$, i.e.
	64	+$m \ge \lceil\log_2 s\rceil$. (The point is not that $\log_2 s$ is a lower bound for all conceivable codes — it is the
	65	+cost of the codeless baseline, in which the message must be self-describing.)
	66	+
	67	+(ii) Lower bound. For any prefix-free code the Kraft inequality holds: $\sum_x 2^{-\ell(x)} \le 1$ (cover2006,
	68	+Thm 5.2.1). Put $q(x) = 2^{-\ell(x)}/Z$ with $Z = \sum_x 2^{-\ell(x)} \le 1$. Then
	69	+$$\sum_x p(x)\ell(x) - H(p) = \sum_x p(x)\log_2\frac{p(x)}{2^{-\ell(x)}} = \sum_x p(x)\log_2\frac{p(x)}{q(x)} - \log_2 Z
	70	+= D_{KL}(p\\|q) - \log_2 Z \ge 0,$$
	71	+because $D_{KL} \ge 0$ (Gibbs' inequality) and $-\log_2 Z \ge 0$. Upper bound. Choose
	72	+$\ell(x) = \lceil -\log_2 p(x)\rceil$. Then $\sum_x 2^{-\ell(x)} \le \sum_x p(x) = 1$, so by the converse Kraft
	73	+theorem (cover2006, Thm 5.2.1) a prefix-free code with these lengths exists, and
	74	+$\sum_x p(x)\ell(x) < \sum_x p(x)(-\log_2 p(x) + 1) = H(p) + 1$. The optimal code does at least as well.
	75	+
	76	+(iii) Subtract (ii) from (i). $\square$
	77	+
	78	+Remark 2.3. Theorem 1 is the source-coding theorem of shannon1948 applied to the plan set. Its content for the
	79	+paper is only this: *a shared code converts a worst-case identification cost into an average-case cost that is lower
	80	+exactly to the extent that tasks are predictable.* The "compression" in the title of the hypothesis refers to this
	81	+step and to nothing else.
	82	+
	83	+---
	84	+
	85	+## 3. Theorem 2 — break-even scale of an institution
	86	+
	87	+DEFINITION 3.1 (institutional cost structure). Adopting a code costs a fixed set-up $F \ge 0$ (designing and
	88	+promulgating the standard), a learning cost $\ell \ge 0$ per agent (each of the $N$ agents must learn the code), and an
	89	+administrative overhead $a \ge 0$ bits per task (verification, record-keeping, enforcement). It saves $G > 0$ bits per
	90	+task (Theorem 1, per recipient; multiply by $k-1$ if desired — the argument is unchanged).
	91	+
	92	+THEOREM 2 (break-even). Over a horizon with $T$ tasks, the institution lowers total communication cost if and only
	93	+if
	94	+$$T\,(G - a) > F + N\ell .$$
	95	+With $T = N\lambda\tau$ over $\tau$ periods, this is equivalent to
	96	+$$\lambda\tau > \frac{F/N + \ell}{G - a}\quad\text{when } G > a,$$
	97	+and the inequality can never hold when $a \ge G$.
	98	+
	99	+Proof. Total cost without the institution is $C_0 = T\lceil\log_2 s\rceil$. With it, $C_1 = F + N\ell + T(L_p^\* + a)$.
	100	+$C_1 < C_0 \iff F + N\ell < T(\lceil\log_2 s\rceil - L_p^\* - a) = T(G - a)$. If $a \ge G$ the right-hand side is
	101	+$\le 0 < F + N\ell$ unless $F = \ell = 0$ and $a = G$, in which case the costs are equal, so strict improvement is
	102	+impossible. If $G > a$, substitute $T = N\lambda\tau$ and divide both sides by $N(G-a) > 0$. $\square$
	103	+
	104	+Corollary 3.2 (scale advantage). The threshold $(F/N + \ell)/(G-a)$ is decreasing in $N$: holding $\lambda,\tau$
	105	+fixed, a larger society crosses the break-even point sooner, because the set-up cost is spread over more agents. It is
	106	+not decreasing in $N$ through $\ell$: per-agent learning costs do not benefit from scale. This is the formal version
	107	+of the observation that institutions "pay off at scale" (johnson1982; garicano2000; bolton1994), and it also makes
	108	+precise when they do not: when overhead $a$ is close to the coding gain $G$, or when horizons $\tau$ are short.
	109	+
	110	+---
	111	+
	112	+## 4. Theorem 3 — mismatched code penalty (institutional rigidity)
	113	+
	114	+THEOREM 3 (cost of a stale code). Let the code be optimal for $p$ in the sense of using lengths
	115	+$\ell_p(x) = \lceil -\log_2 p(x) \rceil$, and let the true task distribution be $q$ with $q \ll p$ (every plan with
	116	+positive $q$-probability has positive $p$-probability). Then the expected length under $q$ satisfies
	117	+$$H(q) + D_{KL}(q\\|p) \le \mathbb{E}_q[\ell_p(X)] < H(q) + D_{KL}(q\\|p) + 1 .$$
	118	+The excess cost of using the $p$-code when the world is $q$, relative to the $q$-optimal code, lies in
	119	+$(D_{KL}(q\\|p) - 1,\ D_{KL}(q\\|p) + 1)$ bits per task.
	120	+
	121	+Proof. Since $-\log_2 p(x) \le \ell_p(x) < -\log_2 p(x) + 1$,
	122	+$$\mathbb{E}_q[\ell_p(X)] \ge \sum_x q(x)\left(-\log_2 p(x)\right) = \sum_x q(x)\log_2\frac{q(x)}{p(x)} + \sum_x q(x)\left(-\log_2 q(x)\right) = D_{KL}(q\\|p) + H(q),$$
	123	+and the same computation with the $+1$ gives the strict upper bound. The $q$-optimal code has expected length in
	124	+$[H(q), H(q)+1)$ by Theorem 1(ii), which yields the stated interval for the excess. (This is the "wrong code" result of
	125	+cover2006, Thm 5.4.3.) $\square$
	126	+
	127	+Corollary 4.1 (fragility). A shock that moves the task distribution from $p$ to $q$ raises the per-task
	128	+communication cost of an existing institution by about $D_{KL}(q\\|p)$ bits. Combined with Theorem 2, the institution
	129	+remains worthwhile after the shock only if $G(p) - D_{KL}(q\\|p) > a$ (up to the $\pm 1$ rounding). An institution can
	130	+therefore lower efficiency after a shock, and the more specialized the code (the smaller $H(p)$ relative to
	131	+$\log_2 s$), the larger the divergence a given shock can produce. This is the precise sense in which standardization
	132	+trades efficiency for rigidity; it echoes tainter1988's declining marginal returns but gives the mechanism a unit (bits
	133	+per task) and a magnitude (a divergence).
	134	+
	135	+---
	136	+
	137	+## 5. Theorem 4 — hierarchical routing bound and root load
	138	+
	139	+DEFINITION 5.1 (hierarchy). A hierarchy is a rooted tree whose $N$ leaves are the agents and whose internal
	140	+nodes are administrative units; every internal node has at most $b \ge 2$ children; the depth (leaf-to-root distance)
	141	+is $d = \lceil \log_b N \rceil$ for a complete $b$-ary tree. A task between two leaves is settled by relay when each
	142	+participant sends a message up to their lowest common ancestor (LCA) and the LCA sends the settlement down.
	143	+
	144	+THEOREM 4.
	145	+(i) Message bound. Any pairwise task can be settled with at most $2d$ messages (up-path plus down-path through the
	146	+LCA), versus up to $N-1$ messages for a flat broadcast in which the initiator must reach every other agent. For $M$
	147	+pairwise tasks the total is at most $2M\lceil\log_b N\rceil$.
	148	+(ii) Root load. In a complete $b$-ary tree with $N = b^d$ leaves, if the two participants of a task are drawn
	149	+uniformly at random among leaves (independently, with replacement), the probability that their LCA is the root is
	150	+$1 - 1/b$. Hence the expected number of tasks that pass through the root is $M(1 - 1/b)$.
	151	+(iii) Federated / modular coordination. If a fraction $\theta \in [0,1]$ of tasks is local — both participants
	152	+belong to the same first-level module (child subtree of the root) by construction of the task process — and the
	153	+remaining $1-\theta$ are drawn uniformly as in (ii), the expected root load is $M(1-\theta)(1 - 1/b)$.
	154	+
	155	+Proof. (i) In a tree there is a unique simple path between two leaves; it passes through their LCA and has length at
	156	+most $2d$ since each leaf is at depth $\le d$ and the LCA is at depth $\ge 0$. Sending one message per edge along
	157	+the two leaf-to-LCA paths settles the task in $\le 2d$ messages. A flat broadcast from one initiator to the other
	158	+$N-1$ agents uses $N-1$ messages. Summing over $M$ tasks gives (i).
	159	+(ii) The root has $b$ children, each the root of a subtree with $b^{d-1}$ leaves. The first participant lies in some
	160	+child subtree $j$; the second participant, drawn uniformly and independently, lies in the same subtree $j$ with
	161	+probability $b^{d-1}/b^d = 1/b$. The LCA is the root exactly when the two participants lie in different child subtrees,
	162	+which has probability $1 - 1/b$. Expectation of a sum of $M$ indicators gives $M(1-1/b)$.
	163	+(iii) Local tasks never reach the root (their LCA lies inside a module); the non-local fraction behaves as in (ii).
	164	+By linearity of expectation the load is $M(1-\theta)(1-1/b)$. $\square$
	165	+
	166	+Remark 5.2. Theorem 4 formalizes the trade-off between centralized administration and federated coordination
	167	+(radner1993; bolton1994): a hierarchy cuts the message cost of coordination from linear to logarithmic in $N$, but
	168	+concentrates a constant fraction $1 - 1/b$ of all non-local tasks on the apex, whose processing rate then bounds $K$
	169	+(Definition 1.3). Modularity reduces apex load by the locality fraction $\theta$ at the price of not reconciling
	170	+across modules. These are the five regimes compared in the pre-registered simulations (no institution, shared norms,
	171	+standardized records, centralized administration, federated), and the simulations measure messages and bytes
	172	+directly rather than assuming the bounds.
	173	+
	174	+---
	175	+
	176	+## 6. Theorem 5 — backlog dynamics under mismatch
	177	+
	178	+DEFINITION 6.1 (backlog). Let $D_t$ be the number of tasks arriving in period $t$ and $K_t$ the number the system
	179	+can settle in period $t$. The backlog of unsettled tasks follows the Lindley recursion
	180	+$$B_{t+1} = \max\{0,\ B_t + D_t - K_t\},\qquad B_0 \ge 0 .$$
	181	+
	182	+THEOREM 5.
	183	+(i) Deterministic case. If $D_t = D$ and $K_t = K$ for all $t$, then: if $D \le K$, $B_t \le \max\{B_0 - t(K-D), 0\}$,
	184	+so the backlog reaches $0$ after at most $\lceil B_0/(K-D) \rceil$ periods when $D < K$ and stays constant at $B_0$
	185	+when $D = K$; if $D > K$, $B_t = B_0 + t(D-K) \to \infty$ linearly.
	186	+(ii) Stochastic case. If $(D_t - K_t)_{t\ge0}$ is i.i.d. with finite mean, then when $\mathbb{E}[D_t] < \mathbb{E}[K_t]$
	187	+the backlog has a proper stationary distribution to which $B_t$ converges from any $B_0$, and when
	188	+$\mathbb{E}[D_t] > \mathbb{E}[K_t]$, $B_t \to \infty$ almost surely. (Standard: the recursion is the single-server
	189	+waiting-time recursion of Lindley (1952), and the stability dichotomy is Loynes' theorem (1962); the proof for (ii) is
	190	+not reproduced here and the result is cited.)
	191	+(iii) Delay under a Markovian approximation. If arrivals are Poisson with rate $\lambda_D$ and settlements are
	192	+exponential with rate $\mu$ (an M/M/1 queue), with $\rho = \lambda_D/\mu < 1$, the expected sojourn time is
	193	+$W = 1/(\mu - \lambda_D) = (1/\mu)/(1-\rho)$, which is increasing and convex in $\rho$ on $[0,1)$ and diverges as
	194	+$\rho \uparrow 1$.
	195	+
	196	+Proof of (i). If $D \le K$, then $B_{t+1} = \max\{0, B_t - (K-D)\}$; by induction
	197	+$B_t \le \max\{0, B_0 - t(K-D)\}$, with equality as long as the argument is nonnegative, which gives the hitting time
	198	+when $D < K$ and constancy when $D = K$. If $D > K$, then $B_t + D - K > 0$ for all $t$ (since $B_t \ge 0$), so the
	199	+maximum is never binding and $B_{t+1} = B_t + (D-K)$, hence $B_t = B_0 + t(D-K)$.
	200	+Proof of (iii). The M/M/1 stationary number in system is geometric with parameter $\rho$, so the mean is
	201	+$\rho/(1-\rho)$ and Little's law gives $W = \rho/((1-\rho)\lambda_D) = 1/(\mu(1-\rho)) = 1/(\mu - \lambda_D)$ (any
	202	+standard queueing reference). The function $f(\rho) = (1-\rho)^{-1}$ has $f' = (1-\rho)^{-2} > 0$ and
	203	+$f'' = 2(1-\rho)^{-3} > 0$ on $[0,1)$, and $f(\rho) \to \infty$ as $\rho \uparrow 1$. $\square$
	204	+
	205	+Scope of Theorem 5. This is a mechanism, not evidence. Polities are not single-server queues: capacity is
	206	+endogenous (institutions are built in response to backlog), arrivals are correlated across periods and agents, and
	207	+"unsettled tasks" are not observed in the historical record. The theorem is used in the paper only to motivate two
	208	+qualitative predictions — that costs rise smoothly and convexly in utilization, and that sustained $\rho > 1$ is
	209	+qualitatively different from $\rho < 1$ — and both are stated as conjectures below.
	210	+
	211	+---
	212	+
	213	+## 7. Conjectures (unproved; each with its falsifier)
	214	+
	215	+CONJECTURE C1 (proxy monotonicity). The pre-registered Seshat and V-Dem proxies are monotone in the model's
	216	+primitives: the scale variables (population, territory, largest settlement, settlement hierarchy) increase with
	217	+$D$ of Definition 1.2, and the four indicator blocks (administration, standardized rules and money, communication
	218	+infrastructure, record-keeping) increase with $K$ of Definition 1.3, so that $P = z(\log D) - z(K)$ is monotone in
	219	+$\log\rho$ up to measurement error. Reason: population and territory raise the number and the dispersion of
	220	+participants; written records, legal codes, money and measurement systems are shared codes (Theorem 1); roads,
	221	+couriers and postal systems raise bits per period; administrative levels implement the routing of Theorem 4.
	222	+Falsifier: the proxies for $D$ and $K$ are so collinear (as Seshat's first principal component suggests; shin2020)
	223	+that $P$ retains almost no variance, or alternative defensible $K$ weightings reverse the sign of $P$ for a
	224	+substantial share of polities — both are reported in `research/measurement.md`.
	225	+
	226	+CONJECTURE C2 (sustained overload raises breakdown probability). Periods with $\rho > 1$ sustained longer than
	227	+the time needed to build capacity raise the probability of institutional breakdown. Mechanism family (not proved):
	228	+unsettled tasks accumulate (Theorem 5), so participants cannot form the common knowledge required for coordinated
	229	+action (chwe2000) and local defection spreads by contagion through the interaction network (morris2000); institutions
	230	+whose codes are stale after a shock lose their cost advantage (Theorem 3) and are abandoned. Falsifier: the
	231	+pre-registered H2 tests find $\beta_P$ not positive, or $P$ does not beat the population-only and capacity-only
	232	+baselines out of sample.
	233	+
	234	+CONJECTURE C3 (smooth, not step). The relation between $P$ and instability is convex and smooth away from
	235	+$\rho \approx 1$ and steepens near it, rather than a discrete threshold. Reason: Theorem 5(iii) gives a smooth,
	236	+convex delay in $\rho$; discontinuities would require an additional tipping mechanism (e.g. a critical mass in a
	237	+contagion model) that the framework does not supply. Falsifier: the pre-registered H4 comparison finds that a hinge
	238	+or step model beats the smooth model out of sample in $\ge 80\%$ of folds, which would favour a tipping mechanism;
	239	+conversely, if neither nonlinear form beats the linear model, C3's "steepening" clause is unsupported.
	240	+
	241	+---
	242	+
	243	+## 8. What is not claimed, and which test speaks to each statement
	244	+
	245	+- No Shannon measurement of history. Theorems 1 and 3 are about codes for a plan set with a known distribution.
	246	+ Seshat indicators (presence of writing, coins, courts, roads) are not message lengths, and no entropy or divergence
	247	+ is computed from historical data anywhere in this project. Shannon quantities appear only in the simulations, where
	248	+ the plan distribution $p$ and the messages are defined by the experimenter and can be measured.
	249	+- No causal claim. The theorems describe an idealized communication system; the empirical tests are associational
	250	+ and pre-registered as such (`research/preregistration.md`, §1 and §5); the identification discussion is in
	251	+ `research/identification.md`.
	252	+- Endogeneity of capacity. The framework takes $K$ as given within a period; historically, $K$ is built in response
	253	+ to $D$ (johnson1982; shin2020). This is why H1 is pre-registered as an association without a direction.
	254	+
	255	+\| Statement \| What can speak to it \| Pre-registered test \|
	256	+\|---\|---\|---\|
	257	+\| Theorem 1 (coding gain) \| Simulated messages and bits with vs. without a shared code \| H3 (regimes b, c vs a) \|
	258	+\| Theorem 2 (break-even) \| Simulated total cost across $N$, $\lambda$, horizon, overhead $a$ \| H3 factorial cells; search for cells where institutions lose \|
	259	+\| Theorem 3 (stale code) \| Simulated shocks to the task distribution \| H3 adverse regimes (shock treatment) \|
	260	+\| Theorem 4 (routing, root load) \| Simulated messages and apex load under centralized vs federated regimes \| H3 (regimes d, e) \|
	261	+\| Theorem 5 (backlog) \| Not testable with the historical data (backlog unobserved); illustrated in simulation only \| none (stated as mechanism) \|
	262	+\| Conjecture C1 \| Measurement diagnostics: collinearity of $D$ and $K$ proxies, sensitivity of $P$ to weights \| `research/measurement.md` \|
	263	+\| Conjecture C2 \| Lagged $P$ and subsequent instability vs baselines, out of sample \| H2 (Seshat Y1–Y3; V-Dem Y4–Y6), H5 \|
	264	+\| Conjecture C3 \| Linear vs hinge vs spline in $P$, out of sample \| H4 \|

added research/measurement.md +80 −0

		`@@ -0,0 +1,80 @@`
	1	+# WP35 — Measurement report (Phase 3)
	2	+
	3	+All numbers come from `results/empirical/measurement_diagnostics.json`, `measurement_diagnostics_modern.json`,
	4	+`measurement_irt_item_parameters.csv`, `power_h2a_*.csv`, produced by `python -m wp35.measurement.build_indices`
	5	+(2026-11-01). No outcome variable was read by that module. Definitions follow `research/preregistration.md` §3.
	6	+
	7	+## 1. Seshat (polity level)
	8	+
	9	+### 1.1 Coordination demand D
	10	+Four logged scale variables (polity population, territory, largest-settlement population; settlement hierarchy levels
	11	+unlogged). First principal component of their pairwise-complete correlation matrix explains 82.1% of variance with
	12	+near-equal loadings (log pop 0.26, log territory 0.24, log capital 0.26, settlement hierarchy 0.24). Score = loading-
	13	+weighted mean of available standardized variables, requiring ≥ 2 of 4; 442 polities have `D_pc1`; 299 have
	14	+population alone (`log_pop`). Range bounds (`value_from`/`value_to`) propagate to `D_pc1_lo/hi`: the median half-width
	15	+is 0.017 SD, the 90th percentile 0.093 SD, so range uncertainty in the scale variables is small relative to
	16	+between-polity variation (the larger uncertainty is which polities have any estimate at all).
	17	+
	18	+### 1.2 Coordination capacity K
	19	+Four blocks of presence indicators, each the share coded present among coded items (≥ 50% coded required):
	20	+administration (5 items incl. rank-scaled administrative levels), standardized rules and money (10), communication
	21	+infrastructure (7), record-keeping (9). `K_mean` = mean of blocks with ≥ 3 blocks: 344 polities. Inter-block
	22	+correlations 0.75–0.79; first PC of the blocks explains 82.4% (loadings 0.25 each), and `K_pc1` correlates 0.999 with
	23	+`K_mean`. Cronbach α of the four blocks by macro-region ranges 0.83 (Europe) to 0.97 (South Asia); mean off-diagonal
	24	+block correlation by era 0.67 (1001–1500) to 0.81 (pre-1000 BCE) — the four capacity blocks hang together
	25	+in every region and era.
	26	+Latent alternative: 2PL IRT (MAP, log-normal prior on discrimination, EM converged in 135 iterations) on the 30 binary
	27	+items (disputed 0.5 → missing; administrative levels excluded) gives `K_irt` for 579 polities (≥ 3 items). It
	28	+correlates 0.96 with `K_mean`. Item difficulties order sensibly (tokens, non-written records, roads, couriers easy;
	29	+paper currency, general postal service, examination system, merit promotion hard; `article` is nearly universal and
	30	+uninformative). Five record-keeping items (written record, script, lists/tables, calendar, practical literature) hit the
	31	+discrimination ceiling (8): they are coded almost jointly in Seshat, i.e. locally dependent; the latent score is
	32	+therefore dominated by the "writing cluster", which is why we keep the transparent block mean as primary. Posterior SD of
	33	+`K_irt`: median 0.34, 90th percentile 0.61 (on a unit-variance scale) — individual-polity capacity scores carry
	34	+substantial uncertainty; this is carried into the econometrics through the IRT variant and the bound runs.
	35	+Cross-region invariance: refitting the IRT separately on the Old World (n = 532) and on the Americas + Southeast
	36	+Asia/Oceania (n = 94) gives item difficulties correlated 0.88 across the two fits (mean absolute difference 1.16
	37	+logits, inflated by boundary estimates on the small New-World sample). Comparability is acceptable for the ordering
	38	+of items but the New-World sample is too small for a strict invariance test; leave-region-out results (H5) are the
	39	+operational safeguard.
	40	+Leave-one-block-out: dropping any block leaves `K` correlated ≥ 0.98 with the full index and `P` correlated ≥ 0.95
	41	+with `P_main`.
	42	+
	43	+### 1.3 Coordination pressure P and the collinearity problem
	44	+`P_main = z(D_pc1) − z(K_mean)` is defined for 317 polities and, in pre-crisis form, for **127 of the 166 CrisisDB
	45	+records**. As anticipated from Turchin et al. (2018) and Shin et al. (2020), scale and capacity are strongly
	46	+collinear: corr(`D_pc1`, `K_mean`) = 0.84 (0.80 with `K_irt`). Consequently `P_main` retains only 16.6% of the
	47	+variance of z(D)+z(K); it correlates 0.39 with D and −0.18 with K, and −0.04 with the Turchin-style single complexity
	48	+axis (`complexity_pc1`, first PC of 57 social-complexity indicators, 49.7% of variance). The mismatch index is thus
	49	+nearly orthogonal to overall complexity, which is exactly what makes the H2 test informative if P carries signal, and
	50	+what makes it fragile: P is a residual, so measurement error in either D or K leaks into it. The pre-registered
	51	+baselines (D-only, K-only, complexity-PC1-only) address this directly.
	52	+
	53	+### 1.4 A-priori power (pre-registration §5.8; random outcomes, actual P distribution, n = 127)
	54	+\| Outcome \| β (per SD of P, logit) \| power \|
	55	+\|---\|---\|---\|
	56	+\| Y1b binary (base rate 0.77) \| 0.4 / 0.5 / 0.6 / 0.8 \| 0.47 / 0.63 / 0.79 / 0.95 \|
	57	+\| Y1 fractional (14 items, base 0.29) \| 0.10 / 0.15 / 0.20 \| 0.49 / 0.81 / 0.96 \|
	58	+Decision (pre-registered rule): the binary severe-outcome test Y1b is underpowered a priori (MDE ≈ 0.6 > 0.5);
	59	+the fractional severity Y1 is adequately powered for effects of ≈ 0.15 log-odds per SD. Y1 remains the primary
	60	+Seshat H2 outcome; Y1b is reported as secondary.
	61	+
	62	+## 2. Modern panel
	63	+`K_modern` = mean of six standardized V-Dem items (≥ 4 coded): rigorous/impartial administration, meritocratic
	64	+appointment, transparent laws, fiscal capacity source, territorial control, EMB capacity; PC1 explains 60.5% with
	65	+loadings 0.13–0.19. `D_modern` = log Maddison population (1820+); `P_modern` defined for 15,697 country-years.
	66	+corr(D, K) = 0.06 (no collinearity problem in the modern panel). Caveat fixed before analysis: corr(`K_modern`,
	67	+`v2x_polyarchy`) = 0.75 — the V-Dem capacity items are strongly associated with electoral democracy even though no
	68	+democracy index enters K; all modern specifications therefore control for polyarchy as pre-registered, and the K-only
	69	+baseline is reported.
	70	+
	71	+## 3. Ethnographic Atlas
	72	+`log_community_size` from EA031 category midpoints (611 societies), `log_population` from continuous EA202 (953),
	73	+`jurisdictional_hierarchy_beyond` = EA033 − 1 (0–4 levels; 1,155). No outcome; used for H1 only, with Glottolog
	74	+family (173 families) and region controls.
	75	+
	76	+## 4. What the measurement model cannot do
	77	+- Seshat K items are presence codes: they measure existence, not reach or quality, of institutions.
	78	+- Most facts are undated, so "pre-crisis" K is mostly whole-polity K (see pre-registration §3.6).
	79	+- No direct indicator of division of labour or dispersion in Polaris; D is scale only.
	80	+- IRT local dependence among writing items; invariance only coarsely checked.

added results/empirical/measurement_diagnostics.json +321 −0

		`@@ -0,0 +1,321 @@`
	1	+{
	2	+ "D_pc1_loadings": {
	3	+ "log_pop": 0.2646,
	4	+ "log_terr": 0.2366,
	5	+ "log_cap": 0.2617,
	6	+ "settl_hier": 0.2371
	7	+ },
	8	+ "D_pc1_var_share": 0.8209,
	9	+ "K_pc1_loadings": {
	10	+ "K_admin": 0.2481,
	11	+ "K_rules_money": 0.2505,
	12	+ "K_comm_infra": 0.2526,
	13	+ "K_records": 0.2488
	14	+ },
	15	+ "K_pc1_var_share": 0.8239,
	16	+ "complexity_pc1_var_share": 0.497,
	17	+ "complexity_pc1_n_vars": 57,
	18	+ "irt_loglik": -3345.5533228338218,
	19	+ "irt_iterations": 135,
	20	+ "n_D_pc1": 442,
	21	+ "n_K_mean": 344,
	22	+ "n_K_irt": 579,
	23	+ "n_P_main": 317,
	24	+ "n_P_main_crisis": 127,
	25	+ "corr_zD_zK": {
	26	+ "D_pc1~K_mean": 0.8384,
	27	+ "D_pc1~K_irt": 0.7998,
	28	+ "log_pop~K_mean": 0.7966,
	29	+ "K_mean~K_irt": 0.9625,
	30	+ "K_mean~K_pc1": 0.9994,
	31	+ "P_main~complexity_pc1": -0.0434,
	32	+ "P_main~D_pc1": 0.3855,
	33	+ "P_main~K_mean": -0.1797
	34	+ },
	35	+ "P_variance_share_of_D_plus_K": 0.1655,
	36	+ "block_correlations": {
	37	+ "K_admin": {
	38	+ "K_admin": 1.0,
	39	+ "K_rules_money": 0.755,
	40	+ "K_comm_infra": 0.766,
	41	+ "K_records": 0.751
	42	+ },
	43	+ "K_rules_money": {
	44	+ "K_admin": 0.755,
	45	+ "K_rules_money": 1.0,
	46	+ "K_comm_infra": 0.789,
	47	+ "K_records": 0.757
	48	+ },
	49	+ "K_comm_infra": {
	50	+ "K_admin": 0.766,
	51	+ "K_rules_money": 0.789,
	52	+ "K_comm_infra": 1.0,
	53	+ "K_records": 0.773
	54	+ },
	55	+ "K_records": {
	56	+ "K_admin": 0.751,
	57	+ "K_rules_money": 0.757,
	58	+ "K_comm_infra": 0.773,
	59	+ "K_records": 1.0
	60	+ }
	61	+ },
	62	+ "block_alpha_by_macro_region": {
	63	+ "Central_Asia": {
	64	+ "n": 17,
	65	+ "alpha": 0.933
	66	+ },
	67	+ "East_Asia": {
	68	+ "n": 35,
	69	+ "alpha": 0.964
	70	+ },
	71	+ "Europe": {
	72	+ "n": 57,
	73	+ "alpha": 0.831
	74	+ },
	75	+ "North_Africa": {
	76	+ "n": 27,
	77	+ "alpha": 0.886
	78	+ },
	79	+ "North_America": {
	80	+ "n": 19,
	81	+ "alpha": 0.936
	82	+ },
	83	+ "South_Asia": {
	84	+ "n": 11,
	85	+ "alpha": 0.968
	86	+ },
	87	+ "Southeast_Asia_Oceania": {
	88	+ "n": 19,
	89	+ "alpha": 0.935
	90	+ },
	91	+ "Southwest_Asia": {
	92	+ "n": 39,
	93	+ "alpha": 0.88
	94	+ }
	95	+ },
	96	+ "block_corr_mean_by_era": {
	97	+ "1-1000": {
	98	+ "n": 46,
	99	+ "mean_offdiag_corr": 0.726
	100	+ },
	101	+ "1000BCE-0": {
	102	+ "n": 35,
	103	+ "mean_offdiag_corr": 0.756
	104	+ },
	105	+ "1001-1500": {
	106	+ "n": 51,
	107	+ "mean_offdiag_corr": 0.666
	108	+ },
	109	+ "1501-1800": {
	110	+ "n": 38,
	111	+ "mean_offdiag_corr": 0.785
	112	+ },
	113	+ "post-1800": {
	114	+ "n": 20,
	115	+ "mean_offdiag_corr": 0.765
	116	+ },
	117	+ "pre-1000BCE": {
	118	+ "n": 43,
	119	+ "mean_offdiag_corr": 0.805
	120	+ }
	121	+ },
	122	+ "measurement_error": {
	123	+ "D_pc1_halfwidth_median": 0.016672931644356626,
	124	+ "D_pc1_halfwidth_q90": 0.09309851215013848,
	125	+ "K_irt_psd_median": 0.34213753777518097,
	126	+ "K_irt_psd_q90": 0.6062046010250374,
	127	+ "sd_K_irt": 1.008347206077521,
	128	+ "P_main_range_median": 0.036226778179016206
	129	+ },
	130	+ "irt_invariance": {
	131	+ "items_compared": 30,
	132	+ "corr_difficulty_old_vs_new": 0.88,
	133	+ "mean_abs_diff_difficulty": 1.162,
	134	+ "items": [
	135	+ "full_time_bureaucrat",
	136	+ "specialized_government_building",
	137	+ "merit_promotion",
	138	+ "examination_system",
	139	+ "formal_legal_code",
	140	+ "judge",
	141	+ "court",
	142	+ "professional_lawyer",
	143	+ "indigenous_coin",
	144	+ "foreign_coin",
	145	+ "paper_currency",
	146	+ "token",
	147	+ "precious_metal",
	148	+ "article",
	149	+ "courier",
	150	+ "postal_station",
	151	+ "general_postal_service",
	152	+ "road",
	153	+ "bridge",
	154	+ "canal",
	155	+ "port",
	156	+ "written_record",
	157	+ "script",
	158	+ "phonetic_alphabetic_writing",
	159	+ "non_phonetic_writing",
	160	+ "lists_tables_and_classification",
	161	+ "calendar",
	162	+ "practical_literature",
	163	+ "nonwritten_record",
	164	+ "mnemonic_device"
	165	+ ],
	166	+ "fits": {
	167	+ "old_world": {
	168	+ "n": 525,
	169	+ "b": [
	170	+ -0.567,
	171	+ -0.598,
	172	+ 0.554,
	173	+ 1.04,
	174	+ -0.659,
	175	+ -0.089,
	176	+ -0.159,
	177	+ 0.605,
	178	+ -0.041,
	179	+ -0.287,
	180	+ 2.037,
	181	+ -1.102,
	182	+ -1.05,
	183	+ -3.99,
	184	+ -0.984,
	185	+ 0.164,
	186	+ 1.134,
	187	+ -1.121,
	188	+ -0.949,
	189	+ -0.328,
	190	+ -0.901,
	191	+ -0.736,
	192	+ -0.845,
	193	+ -0.409,
	194	+ 1.716,
	195	+ -0.673,
	196	+ -0.669,
	197	+ -0.498,
	198	+ -3.397,
	199	+ -1.374
	200	+ ],
	201	+ "a": [
	202	+ 3.44,
	203	+ 3.543,
	204	+ 1.905,
	205	+ 3.328,
	206	+ 2.428,
	207	+ 2.984,
	208	+ 2.47,
	209	+ 2.412,
	210	+ 3.564,
	211	+ 1.758,
	212	+ 1.82,
	213	+ 0.646,
	214	+ 2.036,
	215	+ 0.626,
	216	+ 2.456,
	217	+ 2.305,
	218	+ 1.695,
	219	+ 3.24,
	220	+ 2.977,
	221	+ 2.077,
	222	+ 1.802,
	223	+ 8.0,
	224	+ 8.0,
	225	+ 1.644,
	226	+ 0.822,
	227	+ 8.0,
	228	+ 8.0,
	229	+ 6.448,
	230	+ 1.202,
	231	+ 1.277
	232	+ ]
	233	+ },
	234	+ "new_world_oceania": {
	235	+ "n": 94,
	236	+ "b": [
	237	+ 0.128,
	238	+ 0.085,
	239	+ 2.663,
	240	+ 3.537,
	241	+ 0.354,
	242	+ 0.801,
	243	+ 0.728,
	244	+ 1.073,
	245	+ 1.881,
	246	+ 1.071,
	247	+ 2.593,
	248	+ 1.207,
	249	+ 0.757,
	250	+ -3.221,
	251	+ 0.518,
	252	+ 1.58,
	253	+ 1.786,
	254	+ 0.072,
	255	+ 0.07,
	256	+ 0.565,
	257	+ 0.578,
	258	+ 0.159,
	259	+ 0.004,
	260	+ 0.736,
	261	+ 1.473,
	262	+ 0.642,
	263	+ -0.049,
	264	+ 0.922,
	265	+ -1.575,
	266	+ -0.93
	267	+ ],
	268	+ "a": [
	269	+ 3.407,
	270	+ 1.466,
	271	+ 0.849,
	272	+ 1.541,
	273	+ 3.221,
	274	+ 2.923,
	275	+ 3.622,
	276	+ 3.145,
	277	+ 1.934,
	278	+ 1.517,
	279	+ 1.11,
	280	+ 0.467,
	281	+ 1.44,
	282	+ 0.454,
	283	+ 1.585,
	284	+ 1.687,
	285	+ 1.946,
	286	+ 1.387,
	287	+ 1.46,
	288	+ 1.677,
	289	+ 2.691,
	290	+ 2.461,
	291	+ 2.748,
	292	+ 3.894,
	293	+ 0.629,
	294	+ 2.858,
	295	+ 3.09,
	296	+ 3.732,
	297	+ 0.826,
	298	+ 0.696
	299	+ ]
	300	+ }
	301	+ }
	302	+ },
	303	+ "leave_one_block_out": {
	304	+ "admin": {
	305	+ "corr_K_with_full": 0.9858,
	306	+ "corr_P_with_main": 0.9575
	307	+ },
	308	+ "rules_money": {
	309	+ "corr_K_with_full": 0.9877,
	310	+ "corr_P_with_main": 0.9615
	311	+ },
	312	+ "comm_infra": {
	313	+ "corr_K_with_full": 0.9874,
	314	+ "corr_P_with_main": 0.9657
	315	+ },
	316	+ "records": {
	317	+ "corr_K_with_full": 0.9816,
	318	+ "corr_P_with_main": 0.9474
	319	+ }
	320	+ }
	321	+}
		\ No newline at end of file

added results/empirical/measurement_diagnostics_modern.json +72 −0

		`@@ -0,0 +1,72 @@`
	1	+{
	2	+ "K_modern_pc1_loadings": {
	3	+ "v2clrspct": 0.1859,
	4	+ "v2stcritrecadm": 0.1678,
	5	+ "v2cltrnslw": 0.1906,
	6	+ "v2stfisccap": 0.15,
	7	+ "v2svstterr": 0.1337,
	8	+ "v2elembcap": 0.1721
	9	+ },
	10	+ "K_modern_pc1_var_share": 0.6052,
	11	+ "n_P_modern": 15697,
	12	+ "corr_D_K_modern": 0.0624,
	13	+ "corr_K_modern_polyarchy": 0.7541,
	14	+ "item_corr": {
	15	+ "v2clrspct": {
	16	+ "v2clrspct": 1.0,
	17	+ "v2stcritrecadm": 0.644,
	18	+ "v2cltrnslw": 0.813,
	19	+ "v2stfisccap": 0.473,
	20	+ "v2svstterr": 0.422,
	21	+ "v2elembcap": 0.574
	22	+ },
	23	+ "v2stcritrecadm": {
	24	+ "v2clrspct": 0.644,
	25	+ "v2stcritrecadm": 1.0,
	26	+ "v2cltrnslw": 0.621,
	27	+ "v2stfisccap": 0.438,
	28	+ "v2svstterr": 0.359,
	29	+ "v2elembcap": 0.521
	30	+ },
	31	+ "v2cltrnslw": {
	32	+ "v2clrspct": 0.813,
	33	+ "v2stcritrecadm": 0.621,
	34	+ "v2cltrnslw": 1.0,
	35	+ "v2stfisccap": 0.513,
	36	+ "v2svstterr": 0.424,
	37	+ "v2elembcap": 0.654
	38	+ },
	39	+ "v2stfisccap": {
	40	+ "v2clrspct": 0.473,
	41	+ "v2stcritrecadm": 0.438,
	42	+ "v2cltrnslw": 0.513,
	43	+ "v2stfisccap": 1.0,
	44	+ "v2svstterr": 0.374,
	45	+ "v2elembcap": 0.497
	46	+ },
	47	+ "v2svstterr": {
	48	+ "v2clrspct": 0.422,
	49	+ "v2stcritrecadm": 0.359,
	50	+ "v2cltrnslw": 0.424,
	51	+ "v2stfisccap": 0.374,
	52	+ "v2svstterr": 1.0,
	53	+ "v2elembcap": 0.443
	54	+ },
	55	+ "v2elembcap": {
	56	+ "v2clrspct": 0.574,
	57	+ "v2stcritrecadm": 0.521,
	58	+ "v2cltrnslw": 0.654,
	59	+ "v2stfisccap": 0.497,
	60	+ "v2svstterr": 0.443,
	61	+ "v2elembcap": 1.0
	62	+ }
	63	+ },
	64	+ "corr_K_by_region": {
	65	+ "1": -0.085,
	66	+ "2": 0.259,
	67	+ "3": 0.081,
	68	+ "4": 0.041,
	69	+ "5": -0.045,
	70	+ "6": -0.012
	71	+ }
	72	+}
		\ No newline at end of file

added results/empirical/measurement_irt_item_parameters.csv +31 −0

		`@@ -0,0 +1,31 @@`
	1	+item,a,b,block
	2	+full_time_bureaucrat,3.6527424079806603,-0.41859686853078504,admin
	3	+specialized_government_building,3.066730654084431,-0.5085373165354822,admin
	4	+merit_promotion,1.7802655820815738,0.7021586341368805,admin
	5	+examination_system,3.3869433671617886,1.1375006596571404,admin
	6	+formal_legal_code,2.9525895205555868,-0.4141059941613492,rules_money
	7	+judge,3.552224718836929,0.060069448676890166,rules_money
	8	+court,3.0760342165002372,0.008468451933063378,rules_money
	9	+professional_lawyer,2.72592841935458,0.6535317471947827,rules_money
	10	+indigenous_coin,3.7072105660727472,0.15112395390732436,rules_money
	11	+foreign_coin,1.9763009083036096,-0.08042415203888174,rules_money
	12	+paper_currency,1.4853538625942566,2.2492864209104257,rules_money
	13	+token,0.579698340872784,-0.8294519474758425,rules_money
	14	+precious_metal,2.1558974591242315,-0.6735737580810314,rules_money
	15	+article,0.5286046257791284,-4.364025532653697,rules_money
	16	+courier,2.670690519140207,-0.6724574117368375,comm_infra
	17	+postal_station,2.560978399199805,0.3210643907483727,comm_infra
	18	+general_postal_service,1.9077225274480532,1.1474143346185122,comm_infra
	19	+road,2.8878850939209206,-0.8863284670864775,comm_infra
	20	+bridge,2.7265013330964507,-0.7485140494607205,comm_infra
	21	+canal,2.165331667573359,-0.1798629790826995,comm_infra
	22	+port,2.2796348701360767,-0.5152583130598091,comm_infra
	23	+written_record,8.0,-0.5885046260962389,records
	24	+script,7.974647676175637,-0.6926954931493575,records
	25	+phonetic_alphabetic_writing,2.1001761814441022,-0.1634518168751169,records
	26	+non_phonetic_writing,0.6965894707957191,1.9288626298500053,records
	27	+lists_tables_and_classification,7.643776344260694,-0.4587993518664177,records
	28	+calendar,8.0,-0.5777488445589021,records
	29	+practical_literature,8.0,-0.2566162216282787,records
	30	+nonwritten_record,1.2026821075448688,-2.671075035463131,records
	31	+mnemonic_device,1.038487055773755,-1.398153218842954,records

added results/empirical/power_h2a_binary.csv +8 −0

		`@@ -0,0 +1,8 @@`
	1	+beta,power,n
	2	+0.2,0.1875,127
	3	+0.3,0.31,127
	4	+0.4,0.47,127
	5	+0.5,0.625,127
	6	+0.6,0.785,127
	7	+0.8,0.945,127
	8	+1.0,0.99,127

added results/empirical/power_h2a_fractional.csv +7 −0

		`@@ -0,0 +1,7 @@`
	1	+beta,power,n
	2	+0.1,0.4925,127
	3	+0.15,0.805,127
	4	+0.2,0.9625,127
	5	+0.3,1.0,127
	6	+0.4,1.0,127
	7	+0.5,1.0,127

modified results/logs/execution_log.md +3 −0

		`@@ -17,3 +17,6 @@`
17	17	\| 2026-11-01 03:55 \| 2 \| Glottolog 5.3 languages.csv registered + fetched (language families for Galton controls) \| registry 9 sources \|
18	18	\| 2026-11-01 04:10 \| 2 \| `python -m wp35.harmonization.build_panels` (3 runs: EA033 code parsing fixed via `ord`; EA202 continuous; Y6 redefined as last regime year with irregular end type after inspecting France 1789–1960) \| `data/processed/.parquet`, `results/diagnostics/`, `paper/figures/fig_seshat_temporal_coverage.pdf`, `fig_seshat_missingness.pdf`, `fig_dplace_focal_years.pdf`, `fig_modern_coverage.pdf`, `coverage_summary.json` \|
19	19	\| 2026-11-01 04:15 \| 2 \| Data dictionary + parser tests (`pytest`: 6 passed) \| `data/metadata/data_dictionary.md`, `tests/test_harmonization.py` \|
	20	+\| 2026-11-01 04:40 \| 3 \| Formal framework drafted (fork agent) and reviewed: Theorems 1–5 with proofs, Conjectures C1–C3, mapping to tests; Lindley 1952 / Loynes 1962 verified by Crossref and added \| `research/framework.md` \|
	21	+\| 2026-11-01 04:50 \| 3 \| `python -m wp35.measurement.build_indices` (2 runs; IRT regularized with MAP priors after boundary estimates) \| `data/processed/seshat_indices.parquet`, `seshat_crisis_indices.parquet`, `modern_indices.parquet`; `results/empirical/measurement_`, `power_h2a_.csv`; figures `fig_measurement_DK.pdf`, `fig_irt_invariance.pdf` \|
	22	+\| 2026-11-01 05:00 \| 3 \| Measurement report: corr(D,K)=0.84, P keeps 16.6% of variance; MDE Y1b≈0.6 (underpowered), Y1≈0.15 \| `research/measurement.md` \|

added src/wp35/measurement/build_indices.py +288 −0

		`@@ -0,0 +1,288 @@`
	1	+"""Phase 3 measurement: build D, K, P (and rival covariates) exactly as pre-registered, with alternatives,
	2	+measurement error, cross-region comparability diagnostics and an a-priori power computation.
	3	+Usage: python -m wp35.measurement.build_indices
	4	+No outcome (Y1, Y2, Y3, Y4-6) is read or used anywhere in this module.
	5	+"""
	6	+from __future__ import annotations
	7	+import json
	8	+from pathlib import Path
	9	+import numpy as np
	10	+import pandas as pd
	11	+import matplotlib
	12	+matplotlib.use("Agg")
	13	+import matplotlib.pyplot as plt
	14	+import statsmodels.api as sm
	15	+from wp35.measurement.irt import fit_2pl
	16	+from wp35.harmonization.build_panels import K_BLOCKS, STANDARDS, D_VARS, PROD_INFRA, MILITARY
	17	+
	18	+ROOT = Path(__file__).resolve().parents[3]
	19	+PROC = ROOT / "data" / "processed"; RES = ROOT / "results" / "empirical"; FIG = ROOT / "paper" / "figures"
	20	+RES.mkdir(parents=True, exist_ok=True); FIG.mkdir(parents=True, exist_ok=True)
	21	+SEED = 20261101
	22	+K_ITEMS_MODERN = ["v2clrspct", "v2stcritrecadm", "v2cltrnslw", "v2stfisccap", "v2svstterr", "v2elembcap"]
	23	+BIN_K = [v for b in K_BLOCKS.values() for v in b if v != "administrative_level"]
	24	+
	25	+
	26	+def z(x: pd.Series, ref: pd.Series \| None = None) -> pd.Series:
	27	+ ref = x if ref is None else ref
	28	+ return (x - ref.mean()) / ref.std(ddof=0)
	29	+
	30	+
	31	+def pc1_pairwise(df: pd.DataFrame, min_coded: int = 2):
	32	+ """First principal component of the pairwise-complete correlation matrix; scores = loading-weighted mean of the
	33	+ available standardized variables (requires >= min_coded coded). Loadings signed so that the sum is positive."""
	34	+ Zs = (df - df.mean()) / df.std(ddof=0)
	35	+ C = Zs.corr(min_periods=10).fillna(0).values
	36	+ w, V = np.linalg.eigh(C); v = V[:, -1]; v = v if v.sum() > 0 else -v
	37	+ load = pd.Series(v / np.abs(v).sum(), index=df.columns)
	38	+ avail = Zs.notna()
	39	+ num = (Zs.fillna(0) * load).sum(axis=1); den = (avail * load.abs()).sum(axis=1)
	40	+ score = num / den.where(den > 0)
	41	+ score[avail.sum(axis=1) < min_coded] = np.nan
	42	+ return score, load, float(w[-1] / w.sum())
	43	+
	44	+
	45	+def block_score(wide: pd.DataFrame, items: list[str], min_share: float = 0.5) -> pd.Series:
	46	+ sub = wide.reindex(columns=items).copy()
	47	+ if "administrative_level" in sub.columns:
	48	+ r = sub["administrative_level"].rank(method="average")
	49	+ sub["administrative_level"] = (r - 1) / (r.max() - 1)
	50	+ coded = sub.notna().mean(axis=1)
	51	+ return sub.mean(axis=1).where(coded >= min_share)
	52	+
	53	+
	54	+def build_seshat():
	55	+ wide = pd.read_parquet(PROC / "seshat_polity_wide.parquet")
	56	+ agg = pd.read_parquet(PROC / "seshat_polity_agg.parquet")
	57	+ meta = pd.read_parquet(PROC / "seshat_polity_meta.parquet").set_index("polity_id")
	58	+ logw = agg.pivot(index="polity_id", columns="variable", values="logv_mean")
	59	+ lo = agg.pivot(index="polity_id", columns="variable", values="lo_min"); hi = agg.pivot(index="polity_id", columns="variable", values="hi_max")
	60	+ out = pd.DataFrame(index=wide.index)
	61	+ # ---- D ----
	62	+ Dmat = pd.DataFrame({"log_pop": logw["polity_population"], "log_terr": logw["polity_territory"],
	63	+ "log_cap": logw["population_of_the_largest_settlement"], "settl_hier": wide["settlement_hierarchy"]})
	64	+ out["D_pc1"], D_load, D_var = pc1_pairwise(Dmat)
	65	+ out["log_pop"] = Dmat.log_pop; out["log_terr"] = Dmat.log_terr; out["log_cap"] = Dmat.log_cap; out["settl_hier"] = Dmat.settl_hier
	66	+ # bounds for measurement-error sensitivity (geometric variables only)
	67	+ with np.errstate(divide="ignore"):
	68	+ out["log_pop_lo"] = np.log(lo["polity_population"].where(lo["polity_population"] > 0)); out["log_pop_hi"] = np.log(hi["polity_population"].where(hi["polity_population"] > 0))
	69	+ out["log_terr_lo"] = np.log(lo["polity_territory"].where(lo["polity_territory"] > 0)); out["log_terr_hi"] = np.log(hi["polity_territory"].where(hi["polity_territory"] > 0))
	70	+ # ---- K blocks ----
	71	+ for b, items in K_BLOCKS.items():
	72	+ out[f"K_{b}"] = block_score(wide, items)
	73	+ blocks = out[[f"K_{b}" for b in K_BLOCKS]]
	74	+ out["K_mean"] = blocks.mean(axis=1).where(blocks.notna().sum(axis=1) >= 3)
	75	+ out["K_n_blocks"] = blocks.notna().sum(axis=1)
	76	+ out["K_pc1"], K_load, K_var = pc1_pairwise(blocks, min_coded=3)
	77	+ # IRT on binary items (disputed 0.5 -> missing)
	78	+ Xb = wide.reindex(columns=BIN_K).copy(); Xb = Xb.where(Xb != 0.5)
	79	+ irt = fit_2pl(Xb.values)
	80	+ out["K_irt"] = irt["eap"]; out["K_irt_psd"] = irt["psd"]; out["K_irt_nitems"] = irt["n_items"]
	81	+ item_params = pd.DataFrame({"item": BIN_K, "a": irt["a"], "b": irt["b"], "block": [next(k for k, v in K_BLOCKS.items() if i in v) for i in BIN_K]})
	82	+ item_params.to_csv(RES / "measurement_irt_item_parameters.csv", index=False)
	83	+ # standards sub-index (sparse; latent only)
	84	+ out["K_standards"] = block_score(wide, STANDARDS)
	85	+ # ---- P ----
	86	+ for dname in ["D_pc1", "log_pop"]:
	87	+ for kname in ["K_mean", "K_irt", "K_pc1"]:
	88	+ out[f"P_{dname}_{kname}"] = z(out[dname]) - z(out[kname])
	89	+ out["P_main"] = out["P_D_pc1_K_mean"]
	90	+ # standardization windows (sensitivity)
	91	+ df = out.join(meta[["macro_region", "era"]])
	92	+ out["P_within_era"] = df.groupby("era").D_pc1.transform(z) - df.groupby("era").K_mean.transform(z)
	93	+ out["P_within_region"] = df.groupby("macro_region").D_pc1.transform(z) - df.groupby("macro_region").K_mean.transform(z)
	94	+ out["P_rank"] = out.D_pc1.rank(pct=True) - out.K_mean.rank(pct=True)
	95	+ # bounds
	96	+ Dlo = pd.DataFrame({"log_pop": out.log_pop_lo, "log_terr": out.log_terr_lo, "log_cap": Dmat.log_cap, "settl_hier": Dmat.settl_hier})
	97	+ Dhi = pd.DataFrame({"log_pop": out.log_pop_hi, "log_terr": out.log_terr_hi, "log_cap": Dmat.log_cap, "settl_hier": Dmat.settl_hier})
	98	+ out["D_pc1_lo"] = ((Dlo - Dmat.mean()) / Dmat.std(ddof=0)).fillna(0).mul(D_load).sum(axis=1) / (Dlo.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)
	99	+ out["D_pc1_hi"] = ((Dhi - Dmat.mean()) / Dmat.std(ddof=0)).fillna(0).mul(D_load).sum(axis=1) / (Dhi.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)
	100	+ out.loc[out.D_pc1.isna(), ["D_pc1_lo", "D_pc1_hi"]] = np.nan
	101	+ out["P_main_lo"] = z(out.D_pc1_lo, out.D_pc1) - z(out.K_mean); out["P_main_hi"] = z(out.D_pc1_hi, out.D_pc1) - z(out.K_mean)
	102	+ # drop inferred codes sensitivity: recompute K_mean using only non-inferred facts
	103	+ facts = pd.read_parquet(PROC / "seshat_facts.parquet")
	104	+ ni = facts[(~facts.inferred) & facts.value.notna()].groupby(["polity_id", "variable"]).value.mean().unstack()
	105	+ out["K_mean_noinf"] = pd.concat([block_score(ni, it) for it in K_BLOCKS.values()], axis=1).pipe(lambda b: b.mean(axis=1).where(b.notna().sum(axis=1) >= 3)).reindex(out.index)
	106	+ out["P_main_noinf"] = z(out.D_pc1) - z(out.K_mean_noinf)
	107	+ # ---- rival covariates (excluded from K by design) ----
	108	+ wf = facts[facts.sheet == "Warfare"].groupby(["polity_id", "variable"]).value.mean().unstack()
	109	+ out["X_warfare"] = wf.mean(axis=1).where(wf.notna().mean(axis=1) >= 0.5).reindex(out.index)
	110	+ out["X_prod_infra"] = block_score(wide, PROD_INFRA)
	111	+ mil = wide.reindex(columns=MILITARY).copy(); r = mil["military_level"].rank(); mil["military_level"] = (r - 1) / (r.max() - 1)
	112	+ out["X_military"] = mil.mean(axis=1).where(mil.notna().mean(axis=1) >= 0.5)
	113	+ out["X_moralizing"] = wide.get("moralizing_supernatural_punishment_and_reward")
	114	+ rt = facts[facts.sheet == "Religion"].groupby(["polity_id", "variable"]).value.mean().unstack()
	115	+ out["X_moralizing"] = rt.get("moralizing_supernatural_punishment_and_reward").reindex(out.index) if "moralizing_supernatural_punishment_and_reward" in rt else np.nan
	116	+ # Turchin-style single complexity axis over ALL social-complexity indicators (baseline B3)
	117	+ sc_vars = [v for v in facts[facts.sheet == "Social complexity"].variable.unique() if v in wide.columns]
	118	+ scmat = wide[sc_vars].copy()
	119	+ for v in ["polity_population", "polity_territory", "population_of_the_largest_settlement"]:
	120	+ scmat[v] = logw[v]
	121	+ scmat = scmat.loc[:, scmat.notna().sum() >= 100]
	122	+ out["complexity_pc1"], C_load, C_var = pc1_pairwise(scmat, min_coded=5)
	123	+ out = out.join(meta[["seshat_region", "homeNGA", "macro_region", "era", "start_year", "end_year", "midpoint", "duration", "censored", "centralization", "utm_zone"]])
	124	+ out.index.name = "polity_id"
	125	+ out.reset_index().to_parquet(PROC / "seshat_indices.parquet", index=False)
	126	+ # ---- pre-crisis versions (same loadings/standardization as the full sample) ----
	127	+ pre = pd.read_parquet(PROC / "seshat_precrisis_agg.parquet")
	128	+ pw = pre.pivot(index="crisis_id", columns="variable", values="value_mean"); plog = pre.pivot(index="crisis_id", columns="variable", values="logv_mean")
	129	+ cr = pd.DataFrame(index=pw.index)
	130	+ Dp = pd.DataFrame({"log_pop": plog.get("polity_population"), "log_terr": plog.get("polity_territory"), "log_cap": plog.get("population_of_the_largest_settlement"), "settl_hier": pw.get("settlement_hierarchy")}, index=pw.index)
	131	+ Zp = (Dp - Dmat.mean()) / Dmat.std(ddof=0)
	132	+ cr["D_pc1"] = (Zp.fillna(0) * D_load).sum(axis=1) / (Zp.notna() * D_load.abs()).sum(axis=1).replace(0, np.nan)
	133	+ cr.loc[Zp.notna().sum(axis=1) < 2, "D_pc1"] = np.nan
	134	+ cr["log_pop"] = Dp.log_pop
	135	+ for b, items in K_BLOCKS.items():
	136	+ cr[f"K_{b}"] = block_score(pw, items)
	137	+ bl = cr[[f"K_{b}" for b in K_BLOCKS]]; cr["K_mean"] = bl.mean(axis=1).where(bl.notna().sum(axis=1) >= 3)
	138	+ Xp = pw.reindex(columns=BIN_K).where(lambda d: d != 0.5)
	139	+ # EAP with the full-sample item parameters (no refit): posterior under fixed a, b
	140	+ nodes, weights = np.polynomial.hermite_e.hermegauss(41); weights /= weights.sum()
	141	+ from scipy.special import expit
	142	+ Pq = expit(irt["a"][None, :] * (nodes[:, None] - irt["b"][None, :]))
	143	+ Xv = Xp.values; obs = ~np.isnan(Xv); Xz_ = np.where(obs, Xv, 0)
	144	+ L = Xz_ @ np.log(Pq + 1e-300).T + ((1 - Xz_) * obs) @ np.log(1 - Pq + 1e-300).T + np.log(weights)[None, :]
	145	+ post = np.exp(L - L.max(axis=1, keepdims=True)); post /= post.sum(axis=1, keepdims=True)
	146	+ cr["K_irt"] = np.where(obs.sum(axis=1) >= 3, post @ nodes, np.nan)
	147	+ cr["K_pc1"] = (((bl - blocks.mean()) / blocks.std(ddof=0)).fillna(0) * K_load).sum(axis=1) / (bl.notna() * K_load.abs()).sum(axis=1).replace(0, np.nan)
	148	+ cr.loc[bl.notna().sum(axis=1) < 3, "K_pc1"] = np.nan
	149	+ for dname in ["D_pc1", "log_pop"]:
	150	+ for kname in ["K_mean", "K_irt", "K_pc1"]:
	151	+ cr[f"P_{dname}_{kname}"] = z(cr[dname], out[dname]) - z(cr[kname], out[kname])
	152	+ cr["P_main"] = cr["P_D_pc1_K_mean"]
	153	+ wfp = pre[pre.variable.isin(wf.columns)].pivot(index="crisis_id", columns="variable", values="value_mean")
	154	+ cr["X_warfare"] = wfp.mean(axis=1).where(wfp.notna().mean(axis=1) >= 0.5).reindex(cr.index)
	155	+ cr["X_prod_infra"] = block_score(pw, PROD_INFRA)
	156	+ milp = pw.reindex(columns=MILITARY).copy(); milp["military_level"] = milp["military_level"].rank().pipe(lambda r: (r - 1) / (r.max() - 1)); cr["X_military"] = milp.mean(axis=1).where(milp.notna().mean(axis=1) >= 0.5)
	157	+ scp = pw.reindex(columns=scmat.columns).copy()
	158	+ for v in ["polity_population", "polity_territory", "population_of_the_largest_settlement"]:
	159	+ scp[v] = plog.get(v)
	160	+ Zc = (scp - scmat.mean()) / scmat.std(ddof=0)
	161	+ cr["complexity_pc1"] = (Zc.fillna(0) * C_load).sum(axis=1) / (Zc.notna() * C_load.abs()).sum(axis=1).replace(0, np.nan)
	162	+ cr.loc[Zc.notna().sum(axis=1) < 5, "complexity_pc1"] = np.nan
	163	+ crisis = pd.read_parquet(PROC / "seshat_crisis.parquet").set_index("crisis_id")
	164	+ cr["polity_id"] = crisis.polity_id; cr["crisis_start"] = crisis.crisis_start
	165	+ cr = cr.join(meta[["macro_region", "era", "seshat_region", "homeNGA", "start_year"]], on="polity_id")
	166	+ cr["years_since_start"] = cr.crisis_start - cr.start_year
	167	+ cr.reset_index().to_parquet(PROC / "seshat_crisis_indices.parquet", index=False)
	168	+ # ---- diagnostics: collinearity, comparability, measurement error ----
	169	+ diag = {"D_pc1_loadings": D_load.round(4).to_dict(), "D_pc1_var_share": round(D_var, 4),
	170	+ "K_pc1_loadings": K_load.round(4).to_dict(), "K_pc1_var_share": round(K_var, 4),
	171	+ "complexity_pc1_var_share": round(C_var, 4), "complexity_pc1_n_vars": int(scmat.shape[1]),
	172	+ "irt_loglik": irt["loglik"], "irt_iterations": irt["iterations"],
	173	+ "n_D_pc1": int(out.D_pc1.notna().sum()), "n_K_mean": int(out.K_mean.notna().sum()), "n_K_irt": int(out.K_irt.notna().sum()),
	174	+ "n_P_main": int(out.P_main.notna().sum()), "n_P_main_crisis": int(cr.P_main.notna().sum()),
	175	+ "corr_zD_zK": {"D_pc1~K_mean": round(out.D_pc1.corr(out.K_mean), 4), "D_pc1~K_irt": round(out.D_pc1.corr(out.K_irt), 4),
	176	+ "log_pop~K_mean": round(out.log_pop.corr(out.K_mean), 4), "K_mean~K_irt": round(out.K_mean.corr(out.K_irt), 4),
	177	+ "K_mean~K_pc1": round(out.K_mean.corr(out.K_pc1), 4), "P_main~complexity_pc1": round(out.P_main.corr(out.complexity_pc1), 4),
	178	+ "P_main~D_pc1": round(out.P_main.corr(out.D_pc1), 4), "P_main~K_mean": round(out.P_main.corr(out.K_mean), 4)},
	179	+ "P_variance_share_of_D_plus_K": round(float(out.P_main.var() / (z(out.D_pc1).var() + z(out.K_mean).var())), 4),
	180	+ "block_correlations": blocks.corr().round(3).to_dict(),
	181	+ "block_alpha_by_macro_region": {}, "block_corr_mean_by_era": {},
	182	+ "measurement_error": {"D_pc1_halfwidth_median": float(((out.D_pc1_hi - out.D_pc1_lo) / 2).median()),
	183	+ "D_pc1_halfwidth_q90": float(((out.D_pc1_hi - out.D_pc1_lo) / 2).quantile(0.9)),
	184	+ "K_irt_psd_median": float(out.K_irt_psd.median()), "K_irt_psd_q90": float(out.K_irt_psd.quantile(0.9)),
	185	+ "sd_K_irt": float(out.K_irt.std()), "P_main_range_median": float((out.P_main_hi - out.P_main_lo).abs().median())}}
	186	+ for r, g in df.groupby("macro_region"):
	187	+ b = g[[f"K_{k}" for k in K_BLOCKS]].dropna()
	188	+ if len(b) >= 8:
	189	+ k = b.shape[1]; alpha = k / (k - 1) * (1 - b.var(ddof=1).sum() / b.sum(axis=1).var(ddof=1))
	190	+ diag["block_alpha_by_macro_region"][r] = {"n": int(len(b)), "alpha": round(float(alpha), 3)}
	191	+ for e, g in df.groupby("era"):
	192	+ b = g[[f"K_{k}" for k in K_BLOCKS]]
	193	+ c = b.corr().values; diag["block_corr_mean_by_era"][e] = {"n": int(b.dropna().shape[0]), "mean_offdiag_corr": round(float(c[np.triu_indices(4, 1)].mean()), 3)}
	194	+ # IRT invariance: refit on Old World vs Americas/Oceania and compare difficulties
	195	+ new_world = df.macro_region.isin(["North_America", "South_America", "Southeast_Asia_Oceania"])
	196	+ inv = {}
	197	+ for name, mask in [("old_world", ~new_world), ("new_world_oceania", new_world)]:
	198	+ f = fit_2pl(Xb[mask.values].values); inv[name] = {"n": int(mask.sum()), "b": f["b"].round(3).tolist(), "a": f["a"].round(3).tolist()}
	199	+ bo, bn = np.array(inv["old_world"]["b"]), np.array(inv["new_world_oceania"]["b"])
	200	+ ok = np.isfinite(bo) & np.isfinite(bn) & (np.abs(bo) < 5.9) & (np.abs(bn) < 5.9)
	201	+ diag["irt_invariance"] = {"items_compared": int(ok.sum()), "corr_difficulty_old_vs_new": round(float(np.corrcoef(bo[ok], bn[ok])[0, 1]), 3),
	202	+ "mean_abs_diff_difficulty": round(float(np.abs(bo[ok] - bn[ok]).mean()), 3), "items": BIN_K, "fits": inv}
	203	+ # leave-one-block-out sensitivity of K and P
	204	+ loo = {}
	205	+ for b in K_BLOCKS:
	206	+ others = [f"K_{k}" for k in K_BLOCKS if k != b]; kk = out[others].mean(axis=1).where(out[others].notna().sum(axis=1) >= 2)
	207	+ loo[b] = {"corr_K_with_full": round(float(kk.corr(out.K_mean)), 4), "corr_P_with_main": round(float((z(out.D_pc1) - z(kk)).corr(out.P_main)), 4)}
	208	+ diag["leave_one_block_out"] = loo
	209	+ (RES / "measurement_diagnostics.json").write_text(json.dumps(diag, indent=2))
	210	+ # ---- figures ----
	211	+ fig, axes = plt.subplots(1, 3, figsize=(12, 3.8))
	212	+ axes[0].scatter(out.K_mean, out.D_pc1, s=8, alpha=0.5, c="#444"); axes[0].set_xlabel("K (mean of 4 blocks)"); axes[0].set_ylabel("D (PC1 of log scale variables)"); axes[0].set_title(f"r = {out.D_pc1.corr(out.K_mean):.2f}, n = {int((out.D_pc1.notna() & out.K_mean.notna()).sum())}", fontsize=9)
	213	+ axes[1].scatter(out.K_mean, out.K_irt, s=8, alpha=0.5, c="#444"); axes[1].set_xlabel("K mean"); axes[1].set_ylabel("K IRT (2PL EAP)"); axes[1].set_title(f"r = {out.K_mean.corr(out.K_irt):.2f}", fontsize=9)
	214	+ axes[2].errorbar(item_params.b, range(len(item_params)), xerr=None, fmt="o", ms=3, color="#444"); axes[2].set_yticks(range(len(item_params))); axes[2].set_yticklabels(item_params.item, fontsize=5); axes[2].set_xlabel("IRT difficulty b"); axes[2].set_title("Item difficulties", fontsize=9)
	215	+ fig.tight_layout(); fig.savefig(FIG / "fig_measurement_DK.pdf"); plt.close(fig)
	216	+ fig, ax = plt.subplots(figsize=(4.5, 4.2)); ax.scatter(bo[ok], bn[ok], s=12, c="#444"); lim = [min(bo[ok].min(), bn[ok].min()) - .3, max(bo[ok].max(), bn[ok].max()) + .3]; ax.plot(lim, lim, "k--", lw=.8)
	217	+ ax.set_xlabel("difficulty, Old World fit"); ax.set_ylabel("difficulty, Americas/Oceania fit"); ax.set_title(f"IRT item invariance (r = {diag['irt_invariance']['corr_difficulty_old_vs_new']:.2f})", fontsize=9)
	218	+ fig.tight_layout(); fig.savefig(FIG / "fig_irt_invariance.pdf"); plt.close(fig)
	219	+ return out, cr, diag
	220	+
	221	+
	222	+def build_modern():
	223	+ p = pd.read_parquet(PROC / "modern_panel.parquet")
	224	+ items = p[K_ITEMS_MODERN]
	225	+ Zi = (items - items.mean()) / items.std(ddof=0)
	226	+ p["K_modern"] = Zi.mean(axis=1).where(Zi.notna().sum(axis=1) >= 4)
	227	+ p["K_modern_pc1"], Km_load, Km_var = pc1_pairwise(items, min_coded=4)
	228	+ p["D_modern"] = p.log_pop
	229	+ p["P_modern"] = z(p.D_modern) - z(p.K_modern)
	230	+ p["P_modern_pc1"] = z(p.D_modern) - z(p.K_modern_pc1)
	231	+ p["P_modern_within_year"] = p.groupby("year").D_modern.transform(z) - p.groupby("year").K_modern.transform(z)
	232	+ p["P_modern_area"] = z(((z(p.log_pop) + z(p.log_area)) / 2)) - z(p.K_modern)
	233	+ p.to_parquet(PROC / "modern_indices.parquet", index=False)
	234	+ diag = {"K_modern_pc1_loadings": Km_load.round(4).to_dict(), "K_modern_pc1_var_share": round(Km_var, 4),
	235	+ "n_P_modern": int(p.P_modern.notna().sum()), "corr_D_K_modern": round(float(p.D_modern.corr(p.K_modern)), 4),
	236	+ "corr_K_modern_polyarchy": round(float(p.K_modern.corr(p.v2x_polyarchy)), 4),
	237	+ "item_corr": items.corr().round(3).to_dict(),
	238	+ "corr_K_by_region": {int(r): round(float(g.D_modern.corr(g.K_modern)), 3) for r, g in p.groupby("e_regionpol_6C") if g.P_modern.notna().sum() > 50}}
	239	+ (RES / "measurement_diagnostics_modern.json").write_text(json.dumps(diag, indent=2))
	240	+ return p, diag
	241	+
	242	+
	243	+def power_h2a(cr: pd.DataFrame, n_rep: int = 400):
	244	+ """A-priori MDE for H2-A: random outcomes assigned to the actual P distribution of crisis records (no real Y used)."""
	245	+ rng = np.random.default_rng(SEED)
	246	+ d = cr.dropna(subset=["P_main"]).copy()
	247	+ P = z(d.P_main).values; reg = pd.factorize(d.seshat_region)[0]; n = len(P)
	248	+ base = 0.77 # Y1b positive share read from coverage_summary (122/158); used only to set the intercept
	249	+ res = []
	250	+ for beta in [0.2, 0.3, 0.4, 0.5, 0.6, 0.8, 1.0]:
	251	+ rej = 0
	252	+ for _ in range(n_rep):
	253	+ a0 = np.log(base / (1 - base))
	254	+ y = rng.random(n) < 1 / (1 + np.exp(-(a0 + beta * P)))
	255	+ try:
	256	+ m = sm.Logit(y.astype(float), sm.add_constant(P)).fit(disp=0, cov_type="cluster", cov_kwds={"groups": reg})
	257	+ rej += m.pvalues[1] < 0.05
	258	+ except Exception:
	259	+ pass
	260	+ res.append({"beta": beta, "power": rej / n_rep, "n": n})
	261	+ res = pd.DataFrame(res); res.to_csv(RES / "power_h2a_binary.csv", index=False)
	262	+ # fractional outcome: Binomial(14, logistic(a0 + beta P)) with base severity 0.29
	263	+ res2 = []
	264	+ for beta in [0.1, 0.15, 0.2, 0.3, 0.4, 0.5]:
	265	+ rej = 0
	266	+ for _ in range(n_rep):
	267	+ a0 = np.log(0.29 / 0.71); pr = 1 / (1 + np.exp(-(a0 + beta * P))); y = rng.binomial(14, pr) / 14
	268	+ try:
	269	+ m = sm.GLM(y, sm.add_constant(P), family=sm.families.Binomial()).fit(cov_type="cluster", cov_kwds={"groups": reg})
	270	+ rej += m.pvalues[1] < 0.05
	271	+ except Exception:
	272	+ pass
	273	+ res2.append({"beta": beta, "power": rej / n_rep, "n": n})
	274	+ res2 = pd.DataFrame(res2); res2.to_csv(RES / "power_h2a_fractional.csv", index=False)
	275	+ return res, res2
	276	+
	277	+
	278	+def main():
	279	+ out, cr, diag = build_seshat()
	280	+ p, dm = build_modern()
	281	+ pw1, pw2 = power_h2a(cr)
	282	+ print(json.dumps({k: v for k, v in diag.items() if k in ("D_pc1_loadings", "D_pc1_var_share", "K_pc1_loadings", "n_D_pc1", "n_K_mean", "n_K_irt", "n_P_main", "n_P_main_crisis", "corr_zD_zK", "P_variance_share_of_D_plus_K", "measurement_error", "irt_invariance", "leave_one_block_out", "block_alpha_by_macro_region")}, indent=1, default=str)[:6000])
	283	+ print(json.dumps({k: v for k, v in dm.items() if k != "item_corr"}, indent=1))
	284	+ print("power binary:\n", pw1.to_string(index=False)); print("power fractional:\n", pw2.to_string(index=False))
	285	+
	286	+
	287	+if __name__ == "__main__":
	288	+ main()

added src/wp35/measurement/irt.py +63 −0

		`@@ -0,0 +1,63 @@`
	1	+"""Two-parameter logistic IRT (marginal maximum likelihood, EM with Gauss–Hermite quadrature).
	2	+
	3	+Items are binary (0/1); missing entries contribute nothing to the likelihood (missing-at-random by design, no
	4	+imputation). Returns item discriminations a_j, difficulties b_j, EAP scores and posterior SDs for each unit.
	5	+Deterministic (no randomness). Written for transparency rather than speed (N ≈ 600 × J ≈ 30 is instant).
	6	+"""
	7	+from __future__ import annotations
	8	+import numpy as np
	9	+from scipy.special import expit
	10	+
	11	+
	12	+def fit_2pl(X: np.ndarray, n_quad: int = 41, max_iter: int = 1000, tol: float = 1e-6, min_items: int = 3,
	13	+ prior_sd_log_a: float = 0.5, prior_sd_b: float = 2.0):
	14	+ """MAP variant: log-normal prior on a (mean 0, sd prior_sd_log_a) and normal prior on b (sd prior_sd_b) to keep
	15	+ items that are nearly constant in the sample from running to the boundary (Heywood-type cases)."""
	16	+ X = np.asarray(X, dtype=float)
	17	+ N, J = X.shape
	18	+ obs = ~np.isnan(X)
	19	+ Xz = np.where(obs, X, 0.0)
	20	+ nodes, weights = np.polynomial.hermite_e.hermegauss(n_quad) # standard normal prior
	21	+ weights = weights / weights.sum()
	22	+ a = np.ones(J); b = np.zeros(J)
	23	+ ll_old = -np.inf
	24	+ for it in range(max_iter):
	25	+ # E-step: posterior over nodes for each unit
	26	+ P = expit(a[None, :] * (nodes[:, None] - b[None, :])) # Q x J
	27	+ logP, log1P = np.log(P + 1e-300), np.log(1 - P + 1e-300)
	28	+ # log-lik of unit i at node q = sum_j obs_ij [x_ij logP_qj + (1-x_ij) log(1-P_qj)]
	29	+ L = Xz @ logP.T + ((1 - Xz) * obs) @ log1P.T # N x Q
	30	+ L += np.log(weights)[None, :]
	31	+ m = L.max(axis=1, keepdims=True)
	32	+ post = np.exp(L - m); norm = post.sum(axis=1, keepdims=True); post /= norm
	33	+ ll = float((m.squeeze() + np.log(norm.squeeze())).sum())
	34	+ # M-step: Newton updates per item on expected counts
	35	+ r = post.T @ Xz # Q x J expected correct
	36	+ n = post.T @ obs.astype(float) # Q x J expected tries
	37	+ for j in range(J):
	38	+ for _ in range(10):
	39	+ p = expit(a[j] * (nodes - b[j]))
	40	+ g_a = np.sum((r[:, j] - n[:, j] * p) * (nodes - b[j])) - (np.log(a[j]) / prior_sd_log_a ** 2) / a[j]
	41	+ g_b = np.sum((r[:, j] - n[:, j] * p) * (-a[j])) - b[j] / prior_sd_b ** 2
	42	+ w = n[:, j] * p * (1 - p)
	43	+ H_aa = -np.sum(w * (nodes - b[j]) 2) - (1 - np.log(a[j])) / (prior_sd_log_a 2 * a[j] ** 2)
	44	+ H_bb = -np.sum(w * a[j] 2) - 1 / prior_sd_b 2
	45	+ H_ab = np.sum(w * a[j] * (nodes - b[j]))
	46	+ H = np.array([[H_aa, H_ab], [H_ab, H_bb]]) - 1e-6 * np.eye(2)
	47	+ step = np.linalg.solve(H, np.array([g_a, g_b]))
	48	+ a[j] -= step[0]; b[j] -= step[1]
	49	+ a[j] = float(np.clip(a[j], 0.05, 8.0)); b[j] = float(np.clip(b[j], -6, 6))
	50	+ if np.abs(step).max() < 1e-6:
	51	+ break
	52	+ if abs(ll - ll_old) < tol:
	53	+ break
	54	+ ll_old = ll
	55	+ # final EAP
	56	+ P = expit(a[None, :] * (nodes[:, None] - b[None, :]))
	57	+ L = Xz @ np.log(P + 1e-300).T + ((1 - Xz) * obs) @ np.log(1 - P + 1e-300).T + np.log(weights)[None, :]
	58	+ m = L.max(axis=1, keepdims=True); post = np.exp(L - m); post /= post.sum(axis=1, keepdims=True)
	59	+ eap = post @ nodes
	60	+ psd = np.sqrt(np.maximum(post @ nodes 2 - eap 2, 0))
	61	+ n_items = obs.sum(axis=1)
	62	+ eap = np.where(n_items >= min_items, eap, np.nan); psd = np.where(n_items >= min_items, psd, np.nan)
	63	+ return {"a": a, "b": b, "eap": eap, "psd": psd, "loglik": ll, "iterations": it + 1, "n_items": n_items}
	64

spb/wp35_uqo

WP35 Phase 3: formal framework (5 theorems, 3 conjectures), D/K/P indices with IRT alternative, measurement diagnostics, a-priori power

Theme