Conversation
ALE, MKDA and CBMR all answer a question about convergence; none of them uses the test statistic a study reports alongside its peak. CBES estimates the standardized effect size (Hedges' g) at each voxel from (x, y, z, statistic, N), without imputing effect-size images. The model is a local (kernel-weighted) zero-inflated Tobit likelihood. Reported peaks are a thresholded sample of the effect field, so pooling them naively is a spatial winner's curse; a study that reported nothing contributes the probability that it would have reported nothing. Zero-inflation separates "how many studies have an effect here" (prevalence) from "how big is it in those studies" (g), so that a silent study can be explained as having no effect rather than as evidence for a small common one. Convergence and effect size come out as two parameters of one likelihood rather than as competing answers. Three things mattered more than the choice of estimator, all found by simulation: - Silence is judged over a wider radius than the value-weighting kernel, so a study whose peak landed 6 mm away does not argue against its own effect. - Silent studies are weighted like an average reporting study at that voxel. Weighting them fully, while reporting studies are kernel-discounted, lets silence outvote evidence and pins the estimate ~0.15 below the truth. - Studies that reported nothing at all are recovered from the collection rather than from the coordinates table, which drops them (#294) even though they carry most of the information about prevalence. Also adds create_effect_size_coordinate_studyset, which simulates the reporting process including thresholding and genuine per-study zeros; without a simulator that models thresholding none of this can be validated. Simulated bias at the focus, 30 studies, N in [20, 40]: naive pooling +0.144 at a true g of 0.5, zero-inflated -0.033. Design rationale, the literature it borrows from, measured recovery, and the known gaps are written up in docs/notes/effect_size_cbma.md. Inference is the weakest part and is not yet validated against a reference implementation. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Two validations, and the second one changed the estimator. Against images. The 21 NIDM pain studies have full t maps. Pooling per-study Hedges' g across those images gives a reference; thresholding the same images at p < .001 and keeping only the peaks gives 2,725 foci, 1.2% of voxels, which is all CBES gets. It recovers the reference at rho = 0.84 (r = 0.80). ALE on the identical coordinates gets rho = 0.19 -- not a knock on ALE, which is measuring convergence, but it does say the reported statistics carry nearly all of the effect-size signal and coordinate density almost none of it. The same comparison shows the magnitude is about twice the truth (0.80 against 0.41). Thresholding is not the only selection at work: a reported peak is a local maximum, and its height is inflated on top of the threshold. That is now the dominant bias on real data and the top open item; until it is modelled, g is a relative map. Against a global null. Foci are noise, nothing is there. The parametric standard errors turn out to be unusable: uncorrected p < .05 fired for 40% of voxels (10% with the zero-inflated model) against a nominal 5%, and FDR and Bonferroni rejected somewhere in 100% of null simulations. They were faithfully correcting p-values that were already meaningless -- g/se is not a null-referenced statistic when every peak being pooled was selected for being large. So CBES now takes its uncorrected p-values from a spatial null by default, as ALE and MKDA do: relocate every focus within the mask, refit, and read p off the resulting distribution of |z|. That gives 0.052 and 0.041 against a nominal 0.05, and no corrector rejected anywhere in any null simulation. The stock FDRCorrector and FWECorrector(bonferroni) are valid on top of it, since both are pure functions of the p map; FWECorrector(montecarlo) reuses the maximum-statistic null computed during fit. null_method="parametric" remains available for exploring effect-size maps and warns on construction. Supporting fixes found along the way: - The permutation null now refits under the same selection model as the observed map, via a shared _statistic(). It previously built the null from the naive weighted mean while the observed map came from the Tobit fit. - The spatial kernel is truncated at 1% of its peak. get_ale_kernel keeps every voxel above floating-point zero, so a 10 mm FWHM kernel reached ~26 mm and every voxel in the brain had several "contributing" studies at weights of order 1e-10, which made n_studies meaningless. - Mask images with a trailing singleton volume axis no longer crash. - EM inner loop uses ndtr and converges on prevalence as well as mu, which cuts a whole-brain zero-inflated fit from 39s to 24s. The Monte Carlo null is still the dominant cost and is documented as such. docs/notes/validate_cbes.py reproduces both validations. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
A permutation in CBES is a full refit, so the estimator's speed is its inference budget: a 1000-iteration null cost 6.7 hours on a whole brain. It is now ~1.3 hours, results bit-for-bit unchanged, and cluster-level correction comes out of the same permutations. Where the time went, in order of what it was worth: - Iterate over weighted (study, voxel) pairs rather than the dense study-by-voxel block (4.3x). The EM was evaluating normal CDFs across the full block and then multiplying most of the results by zero: at any voxel a study has either reported nearby or been silent there, and most studies are neither, having reported in the region but outside that voxel's kernel. The arithmetic is identical -- a test pads a problem with studies that say nothing and asserts the answer does not move, and the rewrite was checked against a literal transcription of the dense version to 1e-10. - A scratch bitmap instead of np.unique for deduplicating the voxels a study's coverage spheres reach. - scipy.special.ndtr and an inlined normal density instead of scipy.stats.norm. - Retire converged voxels from the EM working set (a further 15% on null iterations). Voxels converge at very different rates -- 21 iterations pass before even a quarter of them settle -- so iterating the whole block until the slowest one is done wastes most of the work. Cluster-level FWE, on size and on mass, is the correction neuroimaging actually uses and was the obvious gap. Clusters are formed on |z| at the statistic corresponding to voxel_thresh read off the null, not assumed: CBES's z is not standard normal, so a nominal 3.29 is not a p of .001. That threshold is only knowable after permuting, which would mean a second pass; instead a short pilot run fixes it first and fit() records cluster size and mass alongside the voxel-level null in the same refits, so cluster correction costs ~5% more rather than 100% more. Under a global null with 20 simulations, all five available corrections reject in exactly 1 of 20: Bonferroni, FDR, voxel FWE, cluster-size FWE and cluster-mass FWE, against a nominal 0.05. Uncorrected rates are 0.052 at .05 and 0.013 at .01. Voxel-level FWE maps are renamed to ALE's convention (logp_level-voxel, logp_desc-size_level-cluster, logp_desc-mass_level-cluster, each with a signed z companion) so the two estimators can be swapped in a workflow. Also fixes a literal %% leaking into the parametric-null warning. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The validation script fitted at the default cluster_threshold of 0.001 and then asked correct_fwe_montecarlo for a forming threshold of 0.01. Those disagree, so the estimator correctly rebuilt the cluster null rather than reusing the one fit() had already paid for -- doubling the permutations behind every simulation in the sweep. Both now name the same threshold. The reported rates are unaffected: clusters were formed at 0.01 either way. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The zero-inflated Monte Carlo cell was the slowest in the sweep and was committed with an asterisk against the earlier 10-simulation figures. It is now measured: uncorrected rates of 0.042 and 0.009 against nominal 0.05 and 0.01, cluster-level FWE rejecting in 1 of 20, and Bonferroni, FDR and voxel FWE in 2 of 20. That corrects a claim made from the partial run. Not every correction rejects in exactly 1 of 20; three cells are 2 of 20. A true rate of 0.05 produces that 26% of the time, so the run is consistent with correct control but does not demonstrate it, and it does not rule out a mild inflation at the voxel level under the zero-inflated model. Settling that needs a few hundred simulations. Both the memo and the class docstring now say so rather than claiming clean control. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
An image study is the limiting case of a coordinate study: it gives the effect at a voxel with no localization uncertainty and no reporting threshold. It therefore enters the existing likelihood with kernel weight 1 everywhere and never contributes a censoring term, and a collection can mix the two freely (use_images=True). Nothing about the mixture or the EM changes. Mixing exposes an incoherence that coordinates alone hide. Sweeping how many of the 21 NIDM pain studies are supplied as images rather than coordinates moved the pooled effect from 0.848 to 0.447 -- the answer depended on data provenance. The cause is that a reported peak is a local maximum. That bias is larger than previously measured here. The earlier 2x compared a peak to its own study's neighbourhood, which still contains the noise that put the peak there. Against a leave-one-out pooling of the other 20 studies, a reported peak overstates the true local effect roughly fivefold (1.22 against 0.22-0.27). The selection model already removes about half of it. peak_bias=rho treats a reported statistic as measuring g_true/rho and rescales the effect by rho, its variance by rho^2 and the reporting threshold by rho, so the censored likelihood stays coherent and the fit scales exactly linearly in rho. Two calibrations were tried and rejected: the raw peak-to-truth ratio corrects twice, because the selection model got there first; a regression slope of image map on coordinate map is attenuated by voxels where one study peaked and the images say nothing. What works is the ratio of the summaries being compared, fitted on studies supplying both. Split-half on the pain collection gives rho = 0.479 +- 0.056 and takes the held-out coordinate/image ratio from 1.96x to 0.94x, and the provenance spread from 0.401 to 0.045. Validated against a known truth with a random-field generator added to docs/notes/validate_cbes.py, because the existing one cannot test this: it draws a value at the ground-truth location and thresholds it, never selecting a location, so it produces no peak-height bias at all. With real local maxima the correction recovers 0.191 against a true 0.227, from 0.583 uncorrected. Two negative results are recorded in the memo rather than smoothed over. The spatial pattern is not recovered in the sparse-peak regime (r = 0.04 against truth, against 0.50 for images) -- the level is right, the location information is not there. And the protocol fails outright on NiMARE's default NeuroVault studyset, where thresholding leaves 3 of 11 studies with any peak; rho is not a transferable constant either, ranging 0.27-0.62 across reporting thresholds, so it must be calibrated per collection. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The objection to the image-calibrated peak_bias is fair: calibrating it needs images, and a collection with images does not need coordinate-based meta-analysis. The principled alternative is random field theory -- a reported peak is a local maximum of a smooth field, RFT gives its height distribution, so the effect should be deconvolvable from the reported statistic alone. Implemented, validated, and it does not work. The reason is a fact about the data, not the implementation, and it reframes the whole estimator. The machinery is sound. A mixture of RFT null peaks and signal peaks recovers its own parameters on synthetic data (pi0 0.70/0.50/0.90 -> 0.71/0.49/0.90), and the null density matches simulated pure-noise fields to +0.17 z units at fMRI-like smoothness. Modelling a signal peak correctly -- as noncentrality plus peak overshoot, because a local maximum of a smooth field sits above its mean even with no selection, not as N(lambda, 1) -- took the simulated bias from 0.81 to 0.21 against a truth of 0.07. But on real data there is nothing to deconvolve. Reported peaks in the 21 NIDM pain studies average z = 3.639; peaks of pure noise at the same threshold average 3.625. The excess is +0.009, where the effect actually present at those locations would give +1.104. With N ~ 16 no voxel has an appreciable chance of clearing z = 3.29 on signal, so reported peaks sit wherever the noise was largest and their heights are set by the threshold. An ablation confirms it: replacing every reported statistic with a constant, or with a draw from the null peak distribution, moves the correlation with the image-based truth from 0.780 to 0.754 and 0.755. The reported magnitudes are worth 0.026 of correlation; the signal is in where the peaks are and in the sample sizes attached to them. So this ships as a diagnostic rather than a correction. peak_information() reports the excess from the reported statistics and the threshold alone, and CBES warns on fit when it falls below 0.25 z units -- which is the regime where the effect-size scale is not identified from the peaks and no correction computed from them can repair it. The deconvolution itself is kept in docs/notes/peak_height_deconvolution.py, validated, for collections of well-powered studies where the excess would be real. This also explains why the image-calibrated peak_bias works and why it does not transfer: the bias is a near-deterministic function of the reporting threshold and the sample size, so one scalar removes it -- a property of the collection, not of the brain. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Three measurements that answer the practical questions, none of them previously in the repo. How many images calibration needs. Since the fit scales exactly linearly in rho, the relative sampling error of rho is the relative error of the reported effect sizes. Twelve draws per point: 2 images give +-145%, 5 give +-17%, 8 give +-12%, 16 give +-6%. Below five it is not worth applying, which is also why the NeuroVault attempt failed at three. Reported magnitudes track the reporting threshold. Thresholding the same 21 studies at p<.001 and at an FWE-like level gives mean reported |g| of 1.259 and 1.650, against 1.233 and 1.825 predicted by a model containing no effect at all. A paper using FWE correction therefore contributes roughly 1.3x the effect size of one using p<.001, for the same brain. A single scalar rho is not valid across a mixed-threshold literature, and threshold="pooled-min" is wrong there too since it applies the most lenient study's threshold to everyone. What coordinates buy on top of images. Against the 21-study image pooling, 10 images alone reach r=0.893; adding 11 coordinate studies reaches 0.930 corrected, which matches 15 images alone (0.932) -- an exchange rate near two coordinate studies per image. Uncorrected they add +0.01 and push the mean away from the truth. Combined with the ablation, the spatial gain is robust and the magnitude gain is not.
A reported peak is a local maximum that cleared the reporting study's own threshold, so its height is set partly by the effect and partly by (u_k, N_k). `peak_bias` was a single number applied to every study, which assumes the whole literature thresholded alike -- not something a literature search can promise. `peak_bias="per-study"` sets rho_k inversely proportional to `null_peak_mean_g(u_k, N_k)`, the effect size a study would report from a peak of pure noise, normalized so the median reporting study sits at `peak_bias_scale`. That divides out the part of the inflation that varies across a heterogeneous collection, which is the part coordinates alone can identify; the remaining common scale is the one number that still needs images. Thresholds now resolve in three ways beyond the existing keywords: - `threshold="study-min-corrected"` undoes the order statistic. The smallest of m_k peaks above u_k sits above u_k by an amount that grows with m_k, so `infer_threshold_from_minimum` solves for the u_k whose expected minimum is what the study reported. This is for the common case where a paper does not state its threshold and only the lowest significant voxel is available. - `threshold` may name a metadata field, for collections where papers do state it. Studies missing the field take the median of those that have it. - `create_effect_size_coordinate_studyset` accepts a sequence of thresholds and records each study's under `reporting_threshold`, so the heterogeneous case is testable with and without knowing the truth. Internally `_study_thresholds` splits into `_study_cutoffs_z` (per-study thresholds on the z scale), `_peak_bias_factors` (per-study rho) and `_apply_peak_bias` (rescales g, its variance, the censoring threshold and the zero-inflation null variance together), so the rescaling stays coherent. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Truncating the 21 NIDM pain studies to their top 3 peaks moves the apparent threshold from 2.33 to 4.45, and the order-statistic correction only recovers 0.09 of that. Truncation and threshold are not separable from peak heights, so `study-min` and `study-min-corrected` are valid only for papers that tabulate every peak that cleared their threshold; with all peaks kept, both recover the imposed threshold to within 0.01. Say so where a reader will see it, and point at `pooled-min` or the metadata field otherwise. Also update the uninformative-peaks warning, which still said images were the only way to correct the bias. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The previous note cited a top-3-peaks truncation as a reason to prefer `pooled-min`. Papers do not build tables that way: what limits a peak table is one local maximum per cluster, or several separated by 8 mm, plus a cluster extent threshold -- none of which is a global "keep the strongest m". A paper reporting only FWE-surviving clusters has changed its threshold, not truncated, and `study-min` recovers that correctly. The measurement that motivated the caveat used one peak per cluster and recovered every imposed threshold to within 0.01 z, so the realistic case was already passing. Keep the caveat only for a deliberately abridged table, which is genuinely unidentifiable from peak heights, and stop implying it is common. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Between-study differences in the peak-height bias are identified from the coordinates alone, but the one common scale is not: rescaling every coordinate study by the same constant leaves the coordinate-only likelihood unchanged. For a coordinate-only fit that is harmless -- the map is then correct up to a multiplicative constant, which is what a relative effect-size map is. It stops being harmless the moment images join the same fit. Images sit on the true Hedges' g scale, so a mismatched constant makes the two kinds of study disagree about the same voxel -- by a factor of 2.05 on the NIDM pain images -- and the pooled value then depends on how many studies of each kind the collection happens to hold. The previous default of 1.0 was therefore the wrong default whenever images were present, and nothing said so. `peak_bias_scale="auto"` reads the constant off the studies that supplied images: fit the images alone and the coordinates alone, take the ratio over the voxels both cover. The fit is exactly linear in the scale, so a ratio of summaries recovers it, and one pass at 1.0 times the answer is the fit at the answer -- no iteration. Mixing images with the default 1.0 now warns, and asking for "auto" without images falls back to a relative map with a warning rather than failing. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…oes not estimate "Relative map" was doing too much work. Measured, the unknown constant touches the magnitude and nothing else: g, se and g_marginal scale with it to six digits, while z, p and prevalence are unchanged. z divides it out top and bottom, and prevalence is a probability that cancels from the mixture responsibilities. So a coordinate-only fit without images gives valid inference and a valid prevalence map on an absolute scale, with g readable only up to a constant -- about 2.65x high on the NIDM pain collection. That is a limit of the data, not the model: the deconvolution that would identify the scale needs peak heights carrying signal, which underpowered studies do not provide. Also drop the claim that the map is "correct in shape". It recovers the image reference at r = 0.77, and that residual is kernel smearing and the information coordinates lack, not the scale -- two separate limitations that should not be conflated. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Profiling a whole-brain fit (228k voxels, 40 studies, mixed thresholds) puts _censoring_terms at 51% of runtime on its own and 63% counting the _normal_pdf calls it makes. Everything optimized in the previous round -- _accumulate, _coverage_entries, _pool -- is now 1-2s combined, so the EM inner loop is the only thing left worth touching. Only mu moves between EM iterations; sigma and the cutoffs do not. The old code recomputed (c-mu)/sigma, (-c-mu)/sigma, /sigma and /sigma**2 on every one, so five divisions per element per iteration became one multiply and two subtracts with the reciprocals hoisted to the caller. lower is now derived from upper by subtracting 2c/sigma rather than a separate negate-and-divide. Both functions are also written in place. At roughly nine million silent (study, voxel) pairs each temporary is about 72 MB, and the naive expressions allocated some ten of them per call, tens of times per fit and once per Monte Carlo relocation; the pdfs are folded into the limits since neither is needed afterwards. Exact, not approximate: verified identical to 1e-11 relative over 1.2M random cases against the previous implementation. Dropping the lower tail was considered and rejected -- it looks negligible at ndtr(-3.4) ~ 3e-4 but is about a third of the score's numerator. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Loosening the EM convergence test from 1e-5 to 1e-4 is worth 10% and moves no voxel's g by more than 0.001 on a whole-brain fit. Going further is not free: 1e-3 buys another 11% but shifts g by up to 0.031, and 1e-2 is 3.1x at the cost of 0.283. Recording what did not work, so it is not attempted again. float32 for the silent arrays looked like an easy near-2x and is worthless twice over: scipy upcasts float32 inside ndtr, so it is 0.92x -- slower -- and the score suffers catastrophic cancellation in pdf_upper - pdf_lower, with relative error 3e4. The previous commit's hoisting was likewise a wash end to end (395.3s vs 396.9s). Benchmarking says why: ndtr costs 21 ms per million elements against 1.7 ms for exp, so _censoring_terms is two erf evaluations per silent pair per iteration and everything around them is noise. Both are irreducible -- the lower tail carries a third of the score's numerator and cannot be dropped. What remains is not micro-optimization. n_cores is already close to linear on the null and defaults to 1, so it is now documented as the largest speedup a caller can ask for, along with the measurement that a null iteration costs 2.4x the observed fit because relocation scatters foci over more of the brain. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Two commits ago I removed a caveat about truncated peak tables on the grounds that papers do not build tables by keeping the strongest m peaks. That was right, but I generalised it too far: cluster-extent thresholds are just as common and do the same damage, for the same reason. Measured on the 21 NIDM pain studies. Reporting one local maximum per cluster recovers every imposed threshold to within 0.02 z. Adding an extent threshold of k >= 10, 20 or 50 voxels lifts the smallest reported statistic by 0.19, 0.32 and 0.57 z, and the order-statistic correction removes only about a tenth of it -- an extent threshold preferentially drops small clusters, which are the ones with low peaks, so it acts exactly like a stricter height threshold. The consequence is an inferred threshold biased high, silence that looks less surprising than it was, and an under-corrected selection. Recorded on both `infer_threshold_from_minimum` and the `threshold` parameter, with the metadata field named as the way out when the papers state their extent threshold. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The Monte Carlo null refits the whole brain once per relocation, which is where essentially all the runtime goes and why a plain fit costs hours. Profiling shut off the micro-optimisation routes -- ndtr is irreducible and dominates -- so the remaining win had to be algorithmic. Eickhoff et al. (2012) replaced ALE's permutation null with histogram convolution, valid because the ALE statistic is a product over independent studies at each voxel. CBES's statistic is an EM fit, so there is no convolution to do, but the independence that licenses convolution also licenses direct sampling: for study k with m_k foci, each focus lands inside a voxel's kernel support with probability |support|/|mask|, or inside the coverage sphere without reaching the kernel, and a focus that lands is a uniformly chosen one of that study's foci carrying its own g. None of it requires touching the brain. The gain is a change of scaling rather than a constant. Relocation costs n_voxels x n_studies x n_iters; this costs n_draws x n_studies, independent of mask size -- 4e7 study-voxel pairs against 9.1e9 on a whole brain with 40 studies -- and the draws are independent where neighbouring voxels of a relocation share studies. Validated rather than assumed: within 2% of the relocation null's |z| threshold at every p from .05 down to 1e-4, r = 0.999 on the p-values, and rejection fractions matching to three decimals at .01 and .001. Only the uncorrected null comes from here. Familywise error needs the maximum over the brain, a property of exactly the spatial correlation this discards, so correct_fwe_montecarlo still relocates -- the same split ALE uses. One stated approximation: a study contributes a single kernel weight even if two of its foci land near one voxel, negligible at p_kernel ~ 0.006. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…(experimental) Two pieces, at very different levels of confidence. The tail approximation is ready. A permutation p cannot fall below 1/(1+n_iters), so resolving a corrected p of 1e-4 needs ten thousand relocations however uninteresting the other 9999 are. Extreme value theory says exceedances of a high threshold are generalized Pareto whatever the parent, so the tail can be modelled from a few hundred (Winkler et al. 2016). My first version was anticonservative at the far tail -- 1.9e-6 against a 20000-permutation truth of 5e-4 -- for two reasons, both now fixed. The goodness-of-fit test used Cramer-von Mises against a distribution whose parameters came from the same data, where the textbook null does not apply and the test is far too lenient; it now uses a parametric bootstrap. And extrapolating a shape parameter fitted to a few hundred points was unbounded; it now stops two orders below the empirical floor and reports the floor beyond, which is conservative. It also falls back to the empirical tail whenever no fit is accepted, so it can refine a quantized p-value but never manufacture one. The rate-based calibration is not ready, and is opt-in only. It addresses the real finding that the absolute scale *is* identified from coordinates: silence is a probit in g whose slope is pinned by the sample size, so across studies differing in N the pattern of who reported traces a dose-response curve. That identification holds -- fitting g from reporting indicators alone recovers a simulated 0.50 as 0.496 and 0.80 as 0.810, using no peak heights. Crucially the cutoffs and dispersions are left on the true scale here, since rescaling them by rho is precisely what floated two quantities that N had pinned. The estimator built on it does not hold. Pooling one global constant as a ratio of means is dragged down by voxels holding no effect, where the rates correctly say zero but a reported noise peak still converts to a positive g; it returns about 0.5 almost regardless of the truth. So peak_bias_scale="rates" is documented as experimental and "auto" will not select it. Better a labelled experiment than an unvalidated method on the default path. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The rate-based scale failed on real data in the opposite direction from simulation: rho came out 1.0 on the NIDM pain collection against an image-calibrated 0.380, having come out ~0.5 regardless of truth in simulation. One bug, two symptoms. The rate likelihood is a binary regression and only carries curvature where the reporting fraction is intermediate. Where almost nobody reported it is flat near zero, so the fit returns zero whatever the truth -- that is the simulation, dragging rho down. Where almost everybody reported it is monotone with no upper bound, complete separation, so the fit runs to the top of the grid -- that is the pain data, pushing rho up. Pooling a ratio of means across both mixes two regimes in which the estimator is meaningless. So restrict to voxels with a reporting fraction between 0.15 and 0.85, and pool by a slope through the origin rather than a ratio of means, which weights voxels by how much effect they carry instead of letting one with almost none cast an equal vote. When too few voxels qualify it now warns and returns a relative map, rather than silently reporting a clamped 1.0. Still experimental and still not selected by "auto". If this does not land near the independent image calibration, the ratio-of-two-fits shortcut is the wrong estimator and it needs a profile likelihood over the scale -- which uses every voxel weighted by its actual information about rho, instead of a hand-picked subset. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Three attempts at this estimator have now failed the same way. Each fitted two things separately and took a ratio, and a ratio pooled by hand is what broke them: it lets a voxel carrying no effect vote as loudly as one carrying the signal (simulation, rho -> 0.5 whatever the truth), and the rate-only fit it was compared against is meaningless wherever almost nobody or almost everybody reported, the latter being complete separation with an unbounded MLE (the pain collection, rho -> 1.0 against an image-calibrated 0.380). The profile likelihood pools nothing by hand. For each candidate scale the values are rescaled, the model refitted, and the entire log-likelihood evaluated, so a voxel contributes exactly as much as its own information about the scale warrants and an uninformative one contributes nothing. Two details make the comparison across scales legitimate. The censoring term stays on the true scale, cutoffs u_k/sqrt(N_k) and dispersions pinned by N_k, never rescaled -- that is the whole source of identification, and rescaling it too is what made the likelihood flat in the scale to begin with. And rescaling the data carries a Jacobian: p -> c p changes the density by c per observation, so sum(w) log(c) is added back, without which the likelihood is maximised by c -> 0, squeezing the data onto a point and fitting it perfectly. Exposed as peak_bias_scale="mle". Still not selected by "auto", and not described as working, until it is checked against the independent image calibration on real data. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Revalidation against a 40000-permutation reference, three parent distributions, 25 repetitions each. The catastrophic failure is gone -- no more 1.9e-6 against a truth of 5e-4 -- and above the empirical floor the fit is accurate: at p = .05 and .01 it returns 0.85-1.00 of the truth and is essentially never more than twice too small. At and below the floor it is not. For a 500-relocation run, 1/501 = 0.002, and at 0.002 and 0.0005 the fit runs about twice anticonservative in 36-60% of runs on every parent tried. That is precisely the range the method exists to reach, so this implementation does not currently deliver its headline benefit. Rather than ship a claim the measurement contradicts, the fit is now applied only above five times the empirical floor, with the empirical tail kept below. What remains is a smoother, less quantized corrected p-value in the resolvable range -- real but modest -- and no extrapolation past n_iters. An anticonservative familywise p is worse than a quantized one. What would earn the extrapolation back: more exceedances to fit than a few hundred permutations supply, or a better threshold-selection rule than shrinking the tail until a bootstrap goodness-of-fit stops rejecting. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The profile likelihood returns 0.050 on the NIDM pain collection against an image calibration of 0.380 -- and 0.050 is the bottom of its search range, so it did not find a maximum, it ran to the boundary. The reason is structural, and explains the three earlier failures too. The value term of this likelihood is exactly scale-invariant. Rescaling the reported values rescales the re-estimated tau-squared with them, so the log-density gains -log c while the Jacobian adds back +log c; they cancel precisely. The only term that still moves with the scale is the censoring one, and because most study-voxel pairs are silent, that term is maximized by an effect of zero. Nothing opposes it: the reported peak heights carry no scale information in this regime, which is the finding this estimator was built around in the first place. So the information that identifies the scale -- and it does exist, a rate-only fit recovers a simulated 0.50 as 0.496 -- lives in the indicator likelihood prod(1 - P_silent) over the studies that reported, not in the value density the censored likelihood uses for them. Profiling a likelihood that cannot see the reporting rate was never going to recover it, however the pooling was arranged. That is also why every ratio-of-two-fits attempt failed: one side identified the scale, the other was invariant to it, and the pooling broke at whichever extreme dominated. Guard rather than pretend: railing to either end of the range now warns and returns a relative map instead of reporting a grid endpoint as an estimate. The next attempt should treat the reporting rate as a moment condition -- choose the scale whose predicted rate matches the observed one -- which uses the identifying information directly, aggregates over voxels instead of fitting unbounded per-voxel MLEs, and weights nothing by hand. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Profiling the likelihood could not recover the scale because its value term is exactly scale-invariant, leaving only the censoring term, which is maximised by an effect of zero. The identifying information lives in the indicator 1 - P_silent, not the value density, so it is used directly as a moment condition: choose the scale whose predicted reporting rate matches the observed one. Predicted rates rise with the scale and the observed rate is fixed, so the match is unique. Three things follow from why the earlier attempts failed. Match per-study rates, not the overall rate. Prevalence can absorb an overall shift -- more studies having an effect and each effect being larger both raise the rate -- so one number cannot separate them. They separate by sample size: a larger effect raises the rate more for a well-powered study, a larger prevalence raises it uniformly. The per-study rates carry that gradient, which is the same exclusion restriction the identification rests on throughout. Keep the cutoffs and dispersions on the true scale. They are pinned by N, and rescaling them is what floated them and flattened the likelihood to begin with. Aggregate over voxels before matching, so nothing is a per-voxel MLE. That is what removes the complete separation which sent the rate-only fit to its boundary wherever almost every study reported. Evaluated on a 5000-voxel subsample, since these are aggregate rates -- seconds per candidate against the roughly 24 whole-brain fits the profile likelihood needed. Landing on either end of the search range warns and returns a relative map instead of reporting an endpoint as an estimate. Not yet checked against the independent image calibration on real data, and not selected by "auto" until it is. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The censoring term asks whether *this voxel* cleared a study's threshold. The
event that actually happened is whether the study's field reached threshold
anywhere in the region its silence is scored over -- a maximum, not a point. A
20 mm sphere is 4169 voxels at 2 mm, 1.8% of the brain, so the two differ by a
great deal: on the 21 NIDM pain studies the observed coverage rate is 0.677
(range 0.134-0.860), and the pointwise form reaches only 0.440 even at g = 0.8.
It cannot produce the observed rate at any effect size, which is why five
attempts at calibrating the effect-size scale all ran to a search boundary.
Random field theory gives the right quantity directly:
P(reach u) = 1 - exp(-2 E[EC](u - g sqrt(N)))
with E[EC] the resel-weighted Euler characteristic of the region. Both
derivatives are closed form, since the Newton step in the EM needs them.
Three things had to be right, each found by testing rather than reasoning:
Clipping E[EC] at zero to keep it positive breaks the gradient. It flattens the
numerical derivative while leaving the analytic score untouched, so Newton
optimises a function it is not evaluating.
Keeping only the 0- and 3-dimensional terms makes P(silent) non-monotone in the
effect -- a larger effect could make silence likelier -- because the
three-dimensional density rises between z = 1 and z = sqrt(3). All four terms
are kept; the 1- and 2-dimensional ones dominate exactly where the
three-dimensional one misbehaves.
The expansion is invalid below its turning point, counting handles and holes
rather than clusters, and goes negative near z = -1. The argument is clamped
there, giving the smallest non-increasing envelope of the valid branch, with
the derivatives zeroed to match. The turning point is solved for rather than
scanned, since a grid returns the last point still rising and leaves a sliver
of the rising branch inside the envelope.
Not yet wired into the likelihood; that comes next, behind an option, with the
scale calibration retried on top of it.
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
censoring="rft" replaces the pointwise silence probability in the EM with the random-field one: a study is silent when the maximum of its field over the region its silence is scored on fails to clear the threshold, not when a single voxel does. Left opt-in, since this changes the term that makes silence quantitative and so needs its own false-positive validation before it can be a default. The calibration path takes the same model as the fit. Getting that wrong would rebuild the exact bug being fixed -- a regional event compared against a pointwise probability -- so the predicted rates in the moment matching now come from whichever censoring is in force. The active-set compaction carries the regional arrays alongside the pointwise ones, and the censoring closure reads them by name so the EM's shrinking stays correct as they are rebound. One coherence condition, documented rather than hidden: the noncentrality is mu * sqrt(N), so mu has to be on the true effect-size scale. That holds when peak_bias puts the reported values there, and does not when peak_bias is None, where mu still carries the peak inflation and the noncentrality is overstated. The regional term is only meaningful alongside a peak-bias correction. smoothness_fwhm defaults to 8 mm and is a second free constant, which is an honest cost of this change: the predicted rate moves 0.813 / 0.580 / 0.397 at g = 0.1 for 8 / 10 / 12 mm. Unlike rho it is at least bounded to a narrow plausible range and often stated in papers. 18 tests over the changed paths pass; the full suite is still running. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…heir mean The regional term treated a study's noncentrality as exactly mu * sqrt(N). A study's own effect is not mu, it is drawn around it with spread tau, and assuming that away made P(silent) collapse from 0.41 to 4e-5 between g = 0 and g = 0.3 on real data. Silence then becomes near-proof of a null effect, the EM drives mu to zero, the predicted reporting rate stays pinned near its null value and the scale search runs to an end of its range -- which is what it did at every smoothness tried, 8 through 14 mm. The pointwise term this replaced carried that spread all along, through sigma = sqrt(1/N + tau^2). Dropping it was a regression introduced while fixing the regional-versus-pointwise error, not a simplification, and it is the second time in this work that correcting one misspecification quietly broke something the previous code had right. Integrated back with five-node Gauss-Hermite over the study effect. Because the nodes shift with mu while tau does not depend on it, the derivatives are the same quadrature applied to dP/dtheta and d2P/dtheta^2, so the Newton step stays exact -- verified to 8e-11 against finite differences, with monotonicity kept. P(silent) now falls 0.652 to 0.139 across the same range rather than 0.737 to 0.036. Five evaluations per censoring call, on the hot path, and confined to the opt-in regional model. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…recoverable The seventh attempt, quadrature over study effects on top of the regional censoring term, rails like the six before it. Each failure had a real and different cause, each fix was correct on its own terms, and each revealed the next -- which is the pattern of a structural obstacle rather than a sequence of bugs. Written up so the next person does not repeat it. The identification itself is real and is not in doubt: fitting the effect from reporting indicators alone, using no peak heights, recovers a simulated 0.50 as 0.496 and 0.80 as 0.810. Every failure is about extracting that information. Four findings worth keeping. The censored likelihood cannot see the reporting rate, because rescaling the values rescales the re-estimated tau-squared and the Jacobian cancels the rest, leaving its value term exactly scale-invariant. Silence is a regional event, not a pointwise one, and the pointwise form cannot produce the observed coverage rate of 0.677 at any effect size. The regional form then needs a smoothness that coordinates do not carry, and that smoothness decides the answer. And smoothness is not one number: within-brain variation (95/5 ratio 1.99x) exceeds between-study variation (1.6x), an 8x spread in resels against 4x, so a global value mis-predicts reporting in a spatial pattern no single scale can absorb. The conclusion is that converting reporting rates into an absolute effect size requires a local resels-per-voxel map, which is computed from images, and a coordinate-only meta-analysis has none. Only the magnitude is lost: z, the p-values, every corrected map and the prevalence are all invariant to the constant. The three failed options are kept as documented negative results. They warn when they rail and "auto" never selects them. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…SDM's precedent No canonical resels-per-voxel map exists to download: SPM and FSL compute it per analysis from that analysis's residuals, since it reflects acquisition, smoothing and registration rather than anatomy alone. SDM has a close precedent for a different quantity -- correlation templates for grey matter/BOLD, white matter, CSF and fractional anisotropy, built from reference data and freely available (Radua et al., 2014). Those describe how the signal covaries across voxels and shape their imputation kernel; the censoring term here needs how the noise is smoothed, which sets how many effectively independent tests a region holds. Related, not interchangeable. Measured on the 21 pain studies, a template built from 20 explains 20.7% of the spatial smoothness pattern in the study held out (median leave-one-out r 0.456, median pairwise r 0.218). Better than a global constant, which captures none of the pattern by construction, and far from sufficient when resels scale as the inverse cube of FWHM. That figure is a lower bound, attenuated by noise in the held-out study's own local estimate -- which is why the leave-one-out correlation so far exceeds the pairwise one. The better lead is not a template. Radua et al. report their recreation did not depend on FWHM once full anisotropy was used: modelling the spatial structure properly made the nuisance parameter stop mattering rather than requiring a good estimate of it. Recorded as the route worth taking if this is reopened. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…resels Measured on the pain images rather than argued from what the two quantities represent. Across-study signal correlation and local noise smoothness are unrelated: 0.1-0.7% of variance explained at every lag from 2 to 18 mm. The lag had to be checked, because a first pass at 2 mm against ~12 mm smoothness gave a saturated signal measure (median 0.958, 5-95 range 0.878 to 0.985) that could not have correlated with anything. By 6-18 mm the measure has a spread of 0.4 to 0.55 and the relationship is still nil, so the near-zero is real rather than an artefact of the measurement. So the split is clean. SDM's correlation templates do what they were built for and could replace the isotropic Gaussian in the pooling kernel, which is a self-contained improvement worth making on its own merits. They cannot stand in for the noise smoothness the regional censoring term needs. That closes the route which looked most promising after the seven scale attempts. What is left of the lead is narrower: their recreation "did not depend on FWHM when full anisotropy was used", which is a claim about their imputation accuracy, and carrying it over to resel counts is an extrapolation across methods rather than a plan. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Three cuts, each decided by a measurement rather than taste. 3340 lines to 2842. peak_bias_scale loses "rates", "mle" and "rate-match", and with them 462 lines of machinery. Seven attempts at fixing the scale without images all failed, and section 15 records why: the conversion needs a local resels-per-voxel map, which is computed from images. Within-brain smoothness varies more than between-study smoothness, so no global value can stand in. Shipping three broken options is worse than one documented negative result. What remains is a float, or "auto"/"images", which is the route that validated at 0.94x out of sample. threshold collapses "study-min" and "study-min-corrected" into one. Undoing the order statistic is never worse -- it can only lower an over-estimated threshold, and reduces to the raw minimum once a study reports enough peaks -- and the end-to-end comparison gave identical maps. Offering the better and the worse version under different names is not a choice. selection_model loses "tobit". Without a zero component every silence has to be explained as a small shared effect, which lands at -0.16 to -0.34 where the truth is positive. Zero inflation strictly dominates it and there is no regime where the plain version is right. "none" stays as the uncorrected baseline. Two tests reworked rather than deleted: the Tobit comparison becomes a direct test that prevalence and effect separate when half the studies are genuinely null, and the study-min test now checks the inferred cutoff against each study's raw minimum instead of against a removed variant. Still to decide, pending a benchmark: whether censoring needs both options, and whether null_method="parametric" survives now that "approximate" exists -- parametric is known to flag 40% of voxels at p < .05 under a global null and was only ever justified by speed. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…w is The previous se/sd figures were taken while the censoring term was inert, so none of them were measuring the coordinate channel. Re-measured on the field simulator, truth known exactly, 30 replications per arm, stratified by the truth. se/sd is now 1.55 to 3.74 against the 1.1 to 2.1 recorded before -- the interval got worse as the point estimate got better. The excess is located rather than inferred: switching the silence off drops it from 1.78 to 1.29 where the effect is, which is demonstrable only now that the term contributes. So the indicator that corrects the magnitude is the same thing that inflates the error. Coverage is 0.97 to 1.00 in every arm and therefore says nothing: at two images it covers with a half-width of 0.92 of the effect. The docstring now says to read se/sd and the width and never coverage, and gives the widths alongside. Also noted: with selection_model="none" there is no censoring roster, so dof falls back to the Kish count over image weights -- 1 at two images, where a t has a critical value of 12.71. That row's width is a correct statement about two studies rather than a comparable interval. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…iddle ground Two accounts of the -0.039 bias where the effect is largest, both falsified. Flattening the report limb's mu-sensitivity (raising its probability to a power alpha, the probe for "a report is over-credited because a reported peak is a local maximum rather than any exceedance") moved the bias the wrong way and monotonically: -0.061 at alpha 0.1 against -0.039 at alpha 1. And alpha = 1, the principled value, is the rmse optimum -- alpha ~ 2 zeroes the bias and costs 40% of rmse. Extending the report over a small radius, on the reasoning that displacement makes a detecting study name a neighbouring voxel and so contribute nothing, is worse at every radius tried. rmse where the effect is ran 0.070 at the named voxel, 0.113 at 4 mm, 0.122 at 6 mm and 0.126 at 20 mm, and the bias flipped from -0.039 to +0.094 at 4 mm alone: one ring overshoots by more than the original undershoot, because a 4 mm sphere asserts the indicator at seven voxels rather than one. Paired p < 0.0001 in the quiet stratum at every radius. The displacement mechanism is real; the likelihood is simply far more sensitive to asserting |g| >= c where it is false than to discarding an observation. The useful part is that the two probes agree: the shipped configuration is the joint optimum of a two-parameter family rather than an arbitrary choice. The sweep result is recorded on _indicator_entries so it is not re-litigated. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Both are reachable and both carry arithmetic, which is where this module has bitten before. reporting_cutoff_to_g's two-sample branch: the same statistic on the same total N implies about twice the effect, since the groups are half the size. Tested directly rather than inherited, because the last time these designs shared a code path the sample size was reduced by mean instead of sum -- reading [30, 30] as two groups of fifteen and inflating everything downstream by 39%. Also covers the sign, the too-small-study refusals and the design check. reported_minimum_z's t_stat branch: a collection tabulating t must give the same bound as one tabulating the same tail probabilities as z, since the bound is quoted on the z scale. Also covers a table with no statistic column, which is ordinary rather than an error, and one whose column holds nothing usable. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
… was wrong Compression was the earlier design's headline failure -- a reference effect spanning elevenfold came back spanning 1.2-fold -- and it had never been measured for the redesign. The docstring claimed the coordinate channel "corrects the bias rather than the compression". That is measurably wrong and is replaced. Four well-separated foci at true g of 0.2, 0.4, 0.6 and 0.8 inside one map, 8 collections, regressing the estimate on the truth over the signal voxels: g recovers a 4.2-fold range for a true 4-fold (slope 0.729, intercept +0.050), against 3.2-fold for pooling the images alone. So the range is now right and slightly over-spread rather than collapsed. The shape of the improvement matters. g is closer to the truth at the strong foci (0.798 for a true 0.800) and further below at the weak ones (0.188 for a true 0.200, against an inflated 0.242 from the images alone). A focus whose true effect is 0.2 against a cutoff near 0.6 g was reported mostly by luck and should be shrunk -- but that is the whole of the slope of 0.73, so the weak end is where to expect over-correction, not the strong end. This also weakens the "-0.039 bias where the effect is largest" that two earlier experiments chased: a top bin spanning 0.25 to 0.50 of truth averages voxels whose estimate is high with voxels whose estimate is low and reports the mixture as a bias. At the focus itself g is 0.798 for a true 0.800. Slope and intercept are the honest summary, and the docstring now says so for the three-bin tables it still carries. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…I wrote The previous commit read the four foci off their absolute values and concluded "the weak end is where to expect over-correction". The relative errors reverse that. At the weakest focus the images are inflated 21% and g is 6% low, so g is three times more accurate there -- which is precisely what the selection correction exists for: a focus at a true 0.2 against a cutoff near 0.6 g is reported mostly by luck, so pooling only what got reported reads it as far too strong. Mean absolute relative error over the four foci is 5.2% for g against 9.8% for the images alone. The only focus where g is worse is 0.4, 14% low against 12% high. So the slope of 0.729 is not the foci being compressed either. It comes from the blob skirts, where the truth runs 0.05 to 0.2 and both arms are dominated by the floor that reading a map as |g| imposes. The table now carries the relative error beside each value, and says to read the per-focus column rather than the slope. g_marginal is the one map that genuinely is worse at the weak end (-25%, -27%), because prevalence falls toward its floor where few studies reported. Narrowed to "prefer it only when comparing against an image-based reference", which is the case where it is the matching estimand. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The table moved twice because it was measured on 6 and 8 collections, where the standard error on a value of 0.2 to 0.8 is 0.02 to 0.03 -- about 10% relative, the size of the effects being reported. Redone at 24 with the SE shown. What survives: the range. g returns 4.26-fold for a true 4-fold against 3.44-fold for pooling the images alone, so the compression that was the earlier design's headline failure is gone. What does not: "5.2% against 9.8%". Mean absolute relative error over the four foci is 7.5% for g against 6.7% for the images -- a wash, slightly the wrong way. The two arms trade errors focus by focus rather than one dominating, so an averaged number was never the summary. The docstring now says so and gives the per-focus columns. The one clear defect is the 0.4 focus: 18% low, about three standard errors, where the images are 5% high. That is the middle of the detection window (1 of 17 studies reported there), so exactly where the likelihood is most sensitive to the reporting model. Flattening the report limb does not fix it, so it is not the limb's weight. Documented as unexplained and as the first thing to look at. Also recorded on the report-limb comment: the alpha probe has now been re-run on a bed that spans the window, where the first one could not. Its middle stratum had been exactly invariant to alpha -- a dead test rather than a null result, since a single-focus bed's middle stratum is the blob skirt where nothing is reported. alpha = 1 is best at every focus where the limb is live. And g at a focus below the window is exactly invariant to alpha, because nothing is reported and silence is nearly flat in mu there: the coordinate channel contributes nothing usable below detection. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Instrumented the four-foci bed to compare the observed reporting rate against the rate the censored likelihood computes at the true mu, 16 collections: truth observed model ratio 0.2 0.007 0.007 1.00 0.4 0.050 0.089 1.78 0.6 0.373 0.402 1.08 0.8 0.708 0.797 1.13 The likelihood gives a reported voxel P(|g| >= c), but a paper reports it only if it cleared c *and* was a local maximum -- a strictly smaller event. So P(report) is overstated by a factor peaking in the middle of the window and falling to 1.00 below it, where nothing is reported at all. An overstated P(report) against an under-observed count pulls mu down hardest where the factor is largest, which is the whole of the 18% deficit at 0.4. This also retracts a conclusion. The earlier sweep that raised the limb's probability to a power was recorded as falsifying the local-maximum hypothesis; it falsified the probe instead. A power rescales the log-probability uniformly in mu, and the real correction is non-monotone (1.00, 1.78, 1.08, 1.13), so that sweep could not represent the hypothesis it was testing. The radius sweeps stay falsified -- those did test what they claimed. Correcting it needs the survival of a suprathreshold local maximum, which needs a field smoothness this model does not carry; documented as open with the target rates. And states plainly what the same instrument showed about the channels: the indicator channel does not dominate. At the 0.4 focus the count alone implies mu = 0.161 and the images imply about 0.4, and the fit lands at 0.337; below the window the count implies a nonsensical -0.144 and the fit is within 6%. The images carry the magnitude, the indicators move it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Measured against the cleanest reference available: 786 HCP subjects, 480 cut into 16 synthetic studies of 30 with coordinates extracted the way papers produce them, and 306 subjects used to make no coordinate at all giving the truth. So the selection that produces the peaks is statistically independent of the quantity compared against, which no earlier reference for this estimator was. estimate r rank r AUC magnitude two images, pooled +0.845 +0.576 0.973 0.85 g +0.830 +0.575 0.967 0.63 g_marginal +0.785 +0.564 0.943 0.54 Pooling the two images alone wins on every metric, magnitude included. That design has a prevalence of exactly 1 -- every synthetic study draws from the same population -- so there is no between-study absence for the mixture to find and nothing for a selection correction to correct. It does harm anyway: fitted prevalence comes back at 0.664 against a true 1.0, because "failed to clear its threshold" and "has no effect" both explain a silence. The silences push mu down through the censoring term and pi down through the mixture, so g_marginal is shrunk twice and is worst of all. The 21-study pain collection is the other end of the regime: different paradigms and populations, prevalence below 1, and the same correction cuts rmse 23% and bias 47%. Since a user cannot easily tell which end a real collection sits at, the docstring now recommends selection_model="none" for a set of similar studies of the same effect in comparable populations, and says plainly that the estimator offers no diagnostic to settle it -- the quantity that would is prevalence, which is what is wrong when it matters. This is stronger than the existing prevalence caveats, which said only that prevalence is compressed and should be read ordinally. The compression propagates into the magnitude. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The previous commit said the silences push mu down through the censoring term and pi down through the mixture, so g_marginal is shrunk twice. Tested by pinning the prevalence at 1 -- the truth on that bed, and a path the code already supported since _fit_chunk sets pi = 1 whenever the model is not zero-inflated. It makes the magnitude worse, not better: 0.60 of the HCP reference against 0.63, and 0.411 against 0.422 for a true 0.5 on the simulator. Removing the "this study has no effect" escape forces every silence to be explained by a small mu, so mu falls further. The censoring term therefore over-shrinks mu whatever the prevalence does, and a fitted pi below 1 is the model partly absorbing that over-shrinkage rather than compounding it. Which makes the HCP failure the same defect as the overstated P(report) documented above, from another direction: too high a reporting probability against an under-observed count drags mu down, and on a prevalence-1 bed there is no genuine absence for pi to absorb it into. The option is not shipped -- it is never the right choice, and "pin the prevalence" is the obvious-looking fix someone would otherwise reach for. The measurement is recorded so it is not re-tried. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
SDM-PSI takes the same coordinates-plus-maps mixture by design, so it was run on the same HCP studies with the same two supplied as t maps -- parity with what CBES gets. It returns r +0.722, AUC 0.927 and 0.43 of the held-out reference: behind g on every column, and behind the two images alone on every column. That matters for the warning above. Two unrelated methods, given identical data, both do worse with the fourteen coordinate tables than without them. So the finding is a property of thresholded coordinate tables in a prevalence-1 regime rather than a quirk of this estimator's censoring term, which is a stronger claim than CBES alone could support. Also recorded for anyone repeating it: giving SDM the maps closed about half its gap (r +0.563 to +0.722, magnitude 0.21 to 0.43), so most of the apparent difference in the coordinates-only arm was input asymmetry rather than method. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The previous commit said the pain and HCP beds "bracket the regime" -- pain at a prevalence below 1 where the correction helps, HCP at exactly 1 where it hurts. Built the regime as a dial to check it: MOTOR_LH studies that carry the effect mixed with EMOTION_FACES studies that do not, so the true prevalence is designed rather than assumed, with held-out MOTOR subjects fixing both estimands exactly. g recovers 0.68, 0.63, 0.56 and 0.53 of mu at true prevalences of 1.00, 0.75, 0.50 and 0.25. It degrades monotonically as prevalence falls, rather than improving where the correction supposedly has something to correct. So prevalence is not the variable. Nor is the statistic, which was the other candidate: the pain bed reported whole-map rmse while the HCP bed reported a ratio at the truth's top quartile. Measured both on the dial, the arms tie on whole-map rmse at prevalence 1 (0.130 against 0.129) and the images win at 0.50 (0.281 against 0.241) -- so CBES does not win on the pain statistic here either. The beds differ in at least three ways at once (real studies against synthetic, real between-study heterogeneity against essentially none, a study-level reference against a subject-pooled one) and which decides it is unmeasured. The recommendation is now empirical rather than principled: fit it both ways, since selection_model="none" costs a tenth of the runtime, and treat a sharp disagreement as a warning rather than a result. Also replaces the simulator's prevalence-compression figures with the designed ones from this bed -- 0.918, 0.714, 0.590, 0.540 against true 1.00, 0.75, 0.50, 0.25, so it tracks to about 0.5 and then floors near 0.54. A designed prevalence on real subjects is better evidence than a simulator for that claim. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
… not the model Four candidate explanations built as dials on the prevalence bed. Three are rejected: prevalence (g recovers 0.68, 0.63, 0.56, 0.53 of mu at true 1.00, 0.75, 0.50, 0.25 -- degrading monotonically), the choice of statistic (ties or loses on the pain bed's own whole-map rmse), and between-study heterogeneity (rmse 0.126, 0.135, 0.159 at tau 0.0, 0.3, 0.6 against the images' 0.126, 0.123, 0.140 -- no crossing). The fourth moves it, and it is a caution about the pain number rather than an endorsement. The pain reference was an inverse-variance mean of 19 study-level g maps, and Hedges' variance is a function of the observed effect, so a study that drew high gets less weight and such a reference is itself pulled downward. g is pulled downward too. An estimator biased low scores better against a reference biased low. Rebuilding the dial's reference that way -- held-out subjects cut into synthetic reference studies rather than pooled -- takes g from tied (0.126 against 0.126) to winning (0.122 against 0.125), and lifts every ratio (g/mu 0.69 to 0.71, images 0.92 to 0.95). But that is about 2 points of the pain gap's 22, so the direction of the artefact is demonstrated and the magnitude is not. The remaining difference -- real studies against synthetic, with their different scanners, paradigms and sample sizes -- cannot be dialled on this bed, and is now the only untested candidate. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…e route The defect was diagnosed by measurement: the likelihood uses P(|g| >= c) where a paper reports a voxel only if it cleared c *and* was a local maximum, overstating the reporting rate by 1.00, 1.78, 1.08 and 1.13 at true g of 0.2, 0.4, 0.6 and 0.8. Now tested by acting on it. The correction is not a reweighting -- that was tried twice and made every focus worse. It is that both limbs must be complementary probabilities of the same event: with E = P(|g| >= c | mu) and q the chance a study exceeding here names this voxel, a reported pair carries qE and a silent one 1 - qE. Still a proper likelihood, and the arithmetic locates the effect: the report limb's score is (qE)'/(qE) = E'/E, so q cancels, while the silent limb's becomes -qE'/(1 - qE), weakened by roughly q -- exactly the over-shrinkage diagnosed. With q fixed at 0.56, the reciprocal of the 1.78 measured at the 0.4 focus, that focus closes precisely: -11% to +1%. So acting on P(report) moves mu in the predicted direction by the predicted amount, which no amount of reasoning had established. But no constant q helps overall: mean absolute error over the four foci runs 8.9%, 9.4%, 10.3%, 11.0%, 12.2% at q of 1.00, 0.90, 0.80, 0.70, 0.56, because a scalar lifts the weak foci and overshoots the strong ones. The shipped q = 1 is the best constant. The fix therefore needs q(mu), the probability that an exceeding voxel is a local maximum, which is a function of the field's smoothness. Unlike cluster extent, which papers do not report reliably, estimated smoothness is routinely printed by SPM and FSL and stated in methods sections -- so this is a missing input the literature already publishes rather than an unobtainable quantity. Documented as open with that route named. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Two corrections to what the previous two commits shipped. First, the ratio had no uncertainty attached, and it is a ratio of two small rates. With a Poisson interval on the counts totalled over 24 collections: 0.4 gives 1.72 [1.19, 3.03], 0.6 gives 1.00 [0.87, 1.18], 0.8 gives 1.12 [1.01, 1.27], and 0.2 has three reports in total and is unmeasured. So the over-statement at the window is real -- its interval excludes 1 -- but two things were stated with more confidence than the data carry: the value below the window was not measured at all, and q is not demonstrably non-monotone, since the 0.6 and 0.8 intervals overlap. The shape is "well below 1 at the window, at or just past 1 above it", not a peak with two sides. Second, and more important: I recommended computing q(mu) from a reported smoothness via a random-field expected-maxima density. Checked against the measurement, that has the sign backwards. The clump argument gives q ~ 1/clump size ~ u^3 with u = (c - mu)/sigma, which falls as mu rises (u runs +1.10, 0.00, -1.10 across the three foci), while the measured q rises (0.58, 1.00, 0.89). The density describes a zero-mean field and these are signal peaks: at a strong focus the blob's own curvature makes that voxel the local maximum so q approaches 1, while at a marginal focus the noise decides which of several exceeding voxels is named. The governing quantity is the signal's curvature against the noise smoothness, which a coordinate table does not carry and a reported FWHM does not supply. What survives: the defect is real at the window, acting on P(report) closes it, no constant q helps overall, and the correct q needs a quantity this model cannot observe. Harder than a missing metadata field, not easier. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…tables Every warning on this estimator is about *when* to trust the coordinate channel -- the prevalence-1 regime where it does harm, the unexplained disagreement between the pain and HCP beds, the over-shrinkage at the window. None of them can be checked against a given collection. What can be checked is whether the channel is even acting at a voxel, and that was already computed and thrown away: _observed_information accumulates the images' contribution to I_mu_mu and the indicators' contribution separately before summing them, so their ratio costs nothing. At 0 the images carry the estimate alone and the tables changed nothing at that voxel, so none of the coordinate caveats apply there; at 1 the indicators carry it and all of them do. On a 20-study collection with two image donors it runs 0.02 to 1.00, with a median of 0.10 and 0.79 at the focus -- the images carry the estimate almost everywhere and the coordinates take over exactly where studies reported, which is the stratification the design rests on, now readable per voxel instead of only in aggregate. Emitted only under the zero-inflated selection model, since with the selection off the tables are inert and the share is identically zero. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The share was shipped as a stratifier, and the obvious next question was whether it also locates the interval's failure. It does not. Across bands of the share from below 0.05 to above 0.50, se/sd runs 2.93, 3.21, 3.23, 2.97 and 2.28, and the top decile is better than the bottom (2.86 against 3.24): the roughly three-fold over-statement is uniform over the map rather than concentrated in the channel that produces it. This also withdraws an earlier reading. Locating the excess "in the censoring term" rested on comparing fits with selection_model="none", which changes the estimator and the dof fallback together, so it was never a within-fit localisation. The docstring now says so. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…igures Two findings, one of which retracts shipped documentation. The se inflation was never the censoring term. The configuration that settles it is the all-image one: with every study carrying an image the reporting indicator is structurally empty and coordinate_share is identically zero, so the coordinate channel cannot be responsible for anything -- yet that arm measured worst, se/sd of 2.21 at quiet voxels against 1.14 for the same images fitted without the mixture, same bed and same seeds. What it was paying for was a prevalence nothing identified: only the indicator separates "no effect in this study" from "a small effect plus noise", and left free against 20 Gaussian values the mixture returned pi = 0.577 against a true 1.000, whose uncertainty was then profiled out of the information about mu. So pi is now held at 1 wherever no study contributes an indicator, and the Schur complement is skipped there -- the limit cannot be left to the algebra, because with pi clamped just below 1 the cross block tends to a ratio of component densities rather than to zero. The responsibility is forced to exactly 1 so the collapse is an identity rather than a tolerance; the all-image arm now returns the non-mixture fit bit for bit, and a configuration with any indicator at all is untouched. Separately, every se/sd figure on record was computed on sd(|g|). g is a signed inverse-variance mean, so at a quiet voxel the absolute value shrinks the spread to about 0.6 of the real one and puts the mean about 0.8 spreads above zero -- halving the denominator and manufacturing a bias at once. At the foci it is nearly a no-op, which is why the old table's quiet column ran 2.71 to 3.74 while its effect column ran 1.55 to 1.97. Corrected, se/sd is 1.2 to 2.1 and the quiet-stratum bias is ~0.000 rather than +0.045 to +0.114: the estimator is not biased upward where nothing was reported. A textbook all-image inverse-variance arm, which reads 1.00 to 1.24, is now the first row of the arm list so the next such error is caught on sight. Also ruled out as a cause: the reporting rule. The model censors on |g| < c while a paper reports local maxima, and feeding the estimator its own rule instead -- 386 foci per collection against 172 -- moves se/sd from 1.83 to 1.80 and the bias not at all. What remains is a conservative interval where the coordinates act, 1.28 to 2.03, falling monotonically as coordinate_share rises, which is the same mechanism seen from the other side: the share is how much the indicator says about pi. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The previous commit speculated that since the pi/mu ridge runs along curves of roughly constant pi*mu, g_marginal would be the sound thing to bracket. Half of that holds and it is the half that does not help. The product is the steadier estimate -- spread 0.076 at quiet voxels against 0.112 for g -- so the ridge does run where the argument said. But its reported error is worse, 0.349 against 0.204, and se_marginal/sd exceeds se/sd in every stratum but the strongest: 4.58 against 1.83 at quiet voxels, then 2.04, 2.48, 2.43, 1.58. The cause is mechanical. Var(mu*pi) needs the whole 2x2 inverse rather than a Schur complement, so the near-singularity the ridge creates amplifies there instead of cancelling: the determinant is small and the delta method divides by it. So the width needs a better-conditioned variance, not a different estimand, which points at a profile likelihood -- it inverts nothing and needs no degrees of freedom, settling both open questions at once. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Codecov flagged 17 uncovered lines in generate.py. Both new error paths turned out to be covered already; the real gap was prevalence, which appears nowhere in either test file despite being the argument the zero-inflated selection model exists to recover. The two simulators apply the draw in different places -- the field one zeroes the effect before building the map, the point one skips the focus entirely -- so both are checked. The test is behavioural rather than a line-coverage exercise: it asserts that fewer studies report when fewer of them have an effect, and that every study is still present in the collection. A generator that quietly shrank each study's effect instead of zeroing some would pass a coverage check and would make the estimator look good against a truth it was not built for. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…o have The remaining route for the interval was a profile likelihood: it inverts nothing, so it is immune to the near-singularity that degrades both the Schur complement and the delta method, and it needs no dof. Implemented as interval="profile", emitting g_lower and g_upper at roughly the cost of a second fit. It answers the question in a way that disqualifies the question. As pi goes to zero the active component explains nothing and the mixture density tends to the null one at every observation, so the profile log-likelihood has a horizontal asymptote at the null-only value, independent of mu. The interval is therefore bounded if and only if the data reject pi = 0. With a couple of image studies among many tables they almost never do: bounded at 2.5% of voxels at two images, 5.5% at five, 7.2% at ten, and even within 10 mm of the focus only 8.6%, 14.8% and 28.4%. Reparametrising does not help, since pi*mu has the same asymptote along pi -> 0 with mu -> infinity. So se reports a curvature at the point the EM selected, and the likelihood does not support that precision about mu anywhere in this regime. The estimate's stability across replications -- a spread of 0.11 where se says 0.20 -- comes from starting at the pooled image mean and stopping at max_iter, not from the data pinning mu down. That also explains why the se/sd question never resolved: it was asking whether an interval was calibrated for a parameter the data do not bound. P-values are unaffected, coming from the permutation null, and so is the ordering of g. The bound-finding docstring previously asserted the deficit is monotone away from the maximum. It is not, for the reason above, so it now says what the walk actually returns: the edge of the connected component containing the estimate, with the disconnected far branch named rather than ignored. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…aviour The previous commit said the profile interval is bounded "if and only if the data reject pi = 0", which was read off a numeric table and is exactly true only where no study reported at the voxel. Writing each observation's mixture density as pi a_i(mu) + (1 - pi) b_i and taking |mu| to infinity gives the asymptote in closed form: image values and silences both have a_i -> 0 and so drop their mu dependence, but a report has a_i -> 1 and keeps pi in the limit. The docstring now carries that expression, with the pi = 0 likelihood-ratio reading given as the special case it is. The derivation also settles by inspection something that cost an experiment: pi*mu has the same asymptote, approached along pi -> 0 with mu -> infinity, so no reparametrisation of the estimand produces a bounded interval here. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The first result here that came from the algebra before the simulation, and the first that bears on whether to run this estimator on a given collection. A silence constrains (pi, mu) only through P = pi S(mu) + (1 - pi) S(0), one equation in two unknowns, which is the ridge. Studies with different sigma and cutoff supply different equations -- but the cutoff measured in sampling standard deviations is just the reported statistic again, c / sigma about z, so S(0) = 2 Phi(z) - 1 depends on the threshold alone while S(mu) = Phi(z - mu sqrt(n)) - Phi(-z - mu sqrt(n)) depends on both. Varying the threshold moves the two components together through the same tail and leaves the equations nearly collinear; varying the sample size moves the active component through mu sqrt(n) with S(0) exactly fixed, and that contrast is what separates the parameters. Measured with 20 studies in every arm, so none of it is about having more data, using the bounded fraction of the profile interval as the identifiability probe. Against a true prevalence of 1.0: alike studies give 0.422 bounded near signal and a fitted 0.823; spreading the sample size over 12 to 120 gives 0.691 and 0.913; spreading the threshold over 2.3 to 4.5 gives 0.414 and 0.815, which is no change; and spreading both is indistinguishable from spreading the sample size alone. So a magnitude is recoverable from a literature of widely differing sample sizes and not from a literature of uniformly sized studies, however many of them there are. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
…ication The previous commit showed that sample-size spread identifies the magnitude on a bed whose foci run 0.2 to 0.8. The obvious follow-up was that it would also explain why this estimator beats an images-only pool on the pain collection and loses on HCP, since the HCP bed is built from equal-sized studies. It does not, and the reason is worth more than the prediction was. Measured with the subject budget, the study count, the held-out truth and the image studies' own sizes all held fixed, so only the tables' heterogeneity varies and the images-only baseline is identical by construction: spreading the table sizes from a uniform 30 to a range of 9 to 70 moves the recovered magnitude from 0.680 to 0.671 and the fitted prevalence from 0.675 to 0.660, with r and AUC identical to three decimals. The bounded fraction of the profile interval explains why -- it is 0.97 in both arms. HCP's motor signal rejects pi = 0 decisively at the voxels being scored, so mu is already identified and there is no ridge for heterogeneity to break, which is what the derivation predicts for that regime. So identifiability is necessary and not sufficient, and the docstring now says so with the numbers. The useful half is the other inference: with mu identified at 97% of scored voxels, g still recovers 0.68 of the truth where the images recover 0.92, and pi still fits 0.66 against a true 1.000. A flat likelihood cannot be blamed for either, so the over-shrinkage and the under-estimated prevalence are properties of what the censoring term assumes a silence to mean. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Every real-data number here was measured on coordinates re-extracted from each study's map, and that proxy had never been compared against coordinates a human transcribed from a paper. The NIDM pain collection carries both for the same 21 studies, so it can be, and they are not interchangeable: the cluster scheme recovers 23% of pain's published peaks within 8 mm and 42% within 20 mm, and finds 119 peaks where the papers printed 267. A majority of published peaks therefore fall outside even the radius a silence asserts over, so extraction reads more voxels as silent than the literature does. Re-run on the published tables, eight paired splits, with the conventional p < 0.001 as the assumed cut and clamp_threshold lowering it per study to that study's smallest reported statistic: rmse falls 14% against pooling the images alone (p = 0.0024), the mean bias 49% (p = 0.0001), and the +0.137 overestimate at the strongest voxels becomes +0.005, the best-centred top stratum of any arm measured. The ordering is unmoved (rank r -0.004 at p = 0.84, AUC -0.006 at p = 0.30) and the Pearson cost is smaller and no longer significant (0.030, p = 0.075) where on extracted tables it was 0.042 at p = 0.017. So the claim is about the estimator rather than about report_peaks. The docstring now carries both rows, explains why they differ in the direction they do -- more silence, more shrinkage, better rmse on a map dominated by near-zero voxels and a worse-centred top -- and says the published row is the one to quote. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The docstring said a random-field expected-maxima density "gets its sign backwards" at a signal peak and gave the intuition -- the blob's own curvature makes that voxel the local maximum. That is now a derivation, which turns an unexplained failure into a bounded one. Writing the field as Z = m + e with e smooth, stationary and mean-zero, a local maximum needs Z'(0) = 0 and Z''(0) < 0. At a signal peak the first condition reduces to the null one, since m'(0) = 0 there. The second does not: Z''(0) = -kappa + e''(0) with kappa = -m''(0) > 0, so P(Z''(0) < 0) = Phi(kappa / sigma_2). A zero-mean density is the kappa = 0 case, where that is exactly one half, and Phi is strictly increasing in kappa. So the standard density is a lower bound, understating the chance of a maximum by a factor of 2 Phi(kappa / sigma_2) -- one where the mean is flat, tending to two as the peak sharpens. Exactly one of the two conditions defining a maximum carries the signal, and it is the one that density pins at a half. So the limb's error is understood, bounded by two, and not fixable from tables alone: the correction needs a per-study peak sharpness no paper reports. The non-monotone pattern in the measured over-statement is not derived and is not claimed to be. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Profiled a whole-brain fit and found a fifth of it inside gc.collect, which nilearn's safe_get_data runs on every call and masker.transform reaches. At two image studies that is 0.55 of 2.0 seconds; it scales with the image count, so a twenty-image collection spends seconds there. When an image is already on the analysis grid a boolean index gives exactly what transform would, so the round trip buys nothing: 136.7 ms becomes 2.2 ms, a 61.7x saving per read. The fast path is declined when the grids differ, and when the masker carries any signal transformation -- standardisation, detrending, smoothing, a target grid -- because those change the values and skipping them to save time would be a different estimator. Values agree to 1.8e-7, nilearn resampling even onto an identical grid, so the fast path skips an interpolation rather than reproducing it; that is below the float32 the maps are stored in but is not bit-identical. Also corrects a wrong claim in _censoring_terms' docstring. It said the cost was spread over four memory-bound kernels, so what pays is removing a pass. Measured on 600,000 pairs, the two ndtr calls take 10.7 ms and 8.1 ms against 3.0 ms for each density and 0.4 ms for an arithmetic pass: 45% of the time is two normal CDFs, which no surrounding rearrangement reaches. Three were tried -- one reciprocal for two divisions, second-from-first through u phi(u) - l phi(l) = u(phi(u) - phi(l)) + k phi(l), and both -- and gave 1.00x, 1.09x and 1.05x with up to 5e-14 of drift, so none is taken. The lower tail cannot be dropped either: median contribution 2e-5 but maximum 0.30, and above 1% of the score's numerator for 38% of pairs. Under a permutation null this function is 59% of the fit, and the levers that work are n_cores and the existing compaction. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Audited for refactor residue rather than reading for it: an AST pass over every module-level name, class attribute, method, function parameter and dataclass field in the estimator, counting code references rather than textual ones so a docstring mention does not keep something alive. One real finding, and it was a chain of three. _accumulate took a `table` it never touched once the coordinate magnitudes were removed; _pool took one only to forward it there; and _study_cutoffs_z took one it had also stopped reading. All three are private with single call sites, so they are gone. _accumulate's docstring had a sentence explaining why the parameter stayed -- that it "still decides the silence geometry downstream, so the signature does not change" -- which is true of the caller and not of this method, so it went with it. A test calling _pool positionally is updated. Also dropped an unused `small_mask` fixture from one test. Checked and deliberately kept: _preprocess_input, _fit and correct_fwe_montecarlo, which the framework calls rather than this file; self.dataset, which is NiMARE's convention and read in cbma/base.py; _mask_bool_, which is the cache behind a method of nearly the same name; the five attributes read through getattr; se_method="hksj", reachable with selection_model="none"; and the peak_bias/fwhm/stat_column/use_images names, whose only remaining mentions are in a test asserting they are gone. Every field of both dataclasses is read. No unused parameters remain anywhere in the file. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
The CBES docstring had grown to 903 lines, a quarter of the module and 43% of the file in docstrings overall. That is a research log in the wrong place: a user opening help(CBES) needs to know what the maps mean and which to trust, not the autopsy of every route that failed. It is now 205 lines, and the full record -- every table, every retraction, every failed route, and the measurement behind each parameter default -- is reproduced verbatim in notes/cbes-evidence.md in the companion experiments repository, which the docstring points at. What stays is what changes a reader's decisions: how the two channels work, six numbered caveats, and the warnings about when this estimator is worse than pooling the images. Sharpened the guidance on which map to read, which was too blunt. An image-based meta-analysis pools per-study maps, so where a fraction pi of studies carry an effect mu it targets pi*mu -- which makes g_marginal the IBMA-comparable estimand and g a quantity no image-only method estimates. The pain reference was itself a study-pool, so g_marginal's lower rmse there was partly correct estimand matching rather than only shrinkage. What still goes wrong is that the fitted prevalence is biased low, so g_marginal estimates the right quantity with a biased multiplier. Both readings are now stated with the case for each, including that powering a new study wants g, power being defined conditional on the alternative being true. Adds examples/02_meta-analyses/17_plot_cbes.py, in the format of its siblings: simulate a collection with both channels, fit, read the maps, recover the planted effects, choose between the two estimands, plan a study, and the two regimes where this estimator should not be used. Runs end to end. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
Condensing the class docstring left the private helpers still carrying research-log material, found by scanning docstrings for decimals and p-values rather than by eye. Four are trimmed to the decision and its reason, with the numbers moved to notes/cbes-evidence.md: _censoring_terms keeps what the two limbs are and why both are needed, and states that the reported limb is over-stated by a bounded factor and that the function is the estimator's hot spot with 45% of it in two normal CDFs -- the per-kernel timings, the three benchmarked rewrites and the lower-tail measurements are in the evidence file. reporting_cutoff_to_g keeps that the assumed degrees of freedom are load-bearing, with the two figures that make the point, and loses the four-by-four table. _study_cutoffs_z keeps that the threshold cannot be inferred and loses what the retired rule cost. _indicator_entries keeps the asymmetry between a silence over a neighbourhood and a report at one voxel, and the four kinds of omitted pair, which are needed to read the code; it loses the rmse figures behind them. Math derivations stay, since those explain the code rather than justify it. Docstrings are now 30% of the file against 43% before, and the file is 2924 lines against 3636 at the start of this pass. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
It was refused under selection_model="zero-inflated", which reports the censored likelihood's curvature instead, so the only configuration that accepted it was selection_model="none" -- where the coordinate tables are inert and CBES is a plain local random-effects fit of the images. In exactly that case nimare.meta.ibma already offers the same correction as small_sample_correction="knapp-hartung", from an estimator family built for images and not pretending to be a coordinate method. So it was a worse route to a correction that already exists, at the cost of a public parameter, a membership check, a bespoke refusal for the combination that made up most of its surface, 30 lines of implementation and its documentation. Removing it also retires weighted_square, an accumulator in _pool that existed only to feed it -- the kind of cascade worth checking for rather than leaving behind. 94 lines out, 92 tests pass. The derivation is kept in the experiments repository in case it is wanted again. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Closes # .
CBESis a coordinate-based effect-size estimator innimare/meta/cbma/effectsize.py, exported fromnimare.meta.cbma. Where ALE and (M)KDA ask where studies agree, this asks how big the effect is — and it gets there from the coordinates' silence rather than from their reported heights.This description was rewritten from scratch after a redesign. An earlier version of this PR described a different estimator: peak heights pooled as effect sizes, a spatial kernel, a
peak_biascorrection,g_relative/g_absolutemaps, and an arrangement null that shuffled effect sizes among the foci of an analysis. All of that is gone. Anything you may have read aboutpeak_bias,fwhm,stat_column,use_imagesor the arrangement null no longer applies.What it does
A reported peak height is a poor estimate of the effect at its own location: it was selected for being large and it sits wherever noise pushed the maximum. Regressing a held-out truth on the effect size a focus's own table reports gives a slope of 0.08 to 0.18 — one tabulated coordinate explains 5% to 9% of the variance in the effect at its own location. So the heights are not read.
The same selection makes silence informative. If a study reported nothing near a voxel, its effect there did not clear that study's reporting threshold, which is a genuine upper bound. Across a literature those bounds constrain the magnitude.
g/g_varmaps by local inverse-variance weighting with a per-voxel DerSimonian–Lairdtau2. A collection must supply at least one such study; one is enough to run, and a fit without any is refused rather than returned as zeros.coverage_radiusof a voxel is silent there; a study that named the voxel reported at that voxel only. The asymmetry is deliberate — a peak is displaced from the effect, so "someone reported 18 mm away" is not evidence about this voxel — and both limbs are needed, because silences alone leave the model reading the silence fraction against a denominator that excludes every study that reported.mu, reported asg) from how often it is present (pi, reported asprevalence), and reports their product asg_marginal. Nothing is imputed.thresholda float, or a metadata field name for per-study values.clamp_thresholdlowers a study's assumed cut to its own smallest reported statistic when the table contradicts the assumption — a hard upper bound, never below the truth, and bit-identical where it should do nothing.reporting_cutoff_to_gis the only place a reported statistic's scale enters the model.permute-imagesscrambles each image study's values among that study's own voxels and holds the silence pattern fixed, so the censoring term cancels between the observed fit and every permutation.correct_fwe_montecarloadds voxel- and cluster-level FWE with a generalized-Pareto tail.interval="profile"additionally emitsg_lower/g_upperfrom the profile likelihood, which inverts no matrix and needs nodof.analysis_masklets an ROI or partial-coverage study declare which voxels it examined, so its silence elsewhere is not read as evidence.nimare/generate.pygainsn_image_studiesandimage_dironcreate_effect_size_coordinate_studyset, which simulates each study's whole statistic field and reports the local maxima that clear its own threshold — without which there is no way to build a valid CBES input.Maps:
g,se,z,p,logp,tau2,n_studies,n_eff,dof; plusprevalence,g_marginal,se_marginalandcoordinate_shareunder the zero-inflated model; plusg_lower/g_upperunderinterval="profile".Which map to read
An image-based meta-analysis pools per-study effect maps, so where a fraction
piof studies carry an effectmuit targets the population meanpi*mu. That makesg_marginalthe IBMA-comparable estimand, andga quantity no image-only method estimates: the effect among the studies that have it.The choice is not clean, because
piis the weakest part of the fit. On the NIDM pain collection, scored against a pool of held-out studies — itself api*muquantity, sog_marginalis the estimand-matched arm —g_marginalwins on rmse (0.184 against 0.230) but undershoots by 0.291 at the strongest voxels wheregsits at +0.005, because the fitted prevalence is biased low. Sog_marginalestimates the right quantity with a biased multiplier.Read
g_marginalwhen comparing against an IBMA or when studies genuinely differ in whether they carry the effect. Readgwhen the question is how large the effect is where present. To power a new study, useg— power is defined conditional on the alternative being true.What it recovers
On the 21-study NIDM pain collection, split in half so the reference comes from studies the coordinates never touched, eight paired splits, using the coordinates the papers actually printed (the collection's own 267 transcribed peaks):
gg_marginalrmse falls 14% against pooling the images alone (paired p = 0.0024) and the mean bias 49% (p = 0.0001), while the +0.137 overestimate at the strongest voxels becomes +0.005 — the best-centred top stratum measured anywhere in this work. The ordering is unmoved (rank r −0.004 at p = 0.84, AUC −0.006 at p = 0.30) and the Pearson cost is not significant (0.030, p = 0.075). The coordinate channel corrects the level and leaves the pattern alone.
Also run on tables re-extracted from the same studies' maps, which is how every earlier number here was measured. The two are not interchangeable: the cluster scheme recovers only 23% of pain's published peaks within 8 mm and finds 119 peaks where the papers printed 267, so it reads more voxels as silent than the literature does. That arm gives a larger rmse gain (22%) and a worse-centred top (−0.065), consistent with more shrinkage. The published row is the one to quote.
What it does not claim
grecovered 0.63, andprevalenceread 0.66 against a true 1.000. That is what happens when every study really has the effect: there is no absence to find and anypi < 1is error. Which regime a collection is in is not observable from the collection.(pi, mu)through one scalar per distinct (sampling error, cutoff) pair, and because the cutoff in sampling-sd units is the reported statistic again, varying the threshold moves both mixture components together while varying the sample size moves only the active one. At a fixed 20 studies, spreading sample sizes 12 to 120 lifts the fraction of voxels where the likelihood boundsmufrom 0.42 to 0.69 and the fitted prevalence from 0.82 to 0.91; spreading thresholds 2.3 to 4.5 changes neither.seis conservative and the likelihood often does not boundmuat all.se/sdruns 1.2 to 2.2 against a known truth, most of it the cost of profiling outpi. The profile interval is unbounded wherever the data do not rejectpi = 0, which at two image studies is 97% of voxels — the honest answer, not a search failure. Read the width, never coverage: every configuration measured covers 0.95 to 1.00, including those admitting almost any magnitude.prevalenceis ordinal, not a fraction, and not comparable between voxels of one map.gvaries 38% over plausible values.max_iterdecides which point on it is returned, reproducibly and no more.Error rates: re-measurement in progress
Flagging this rather than burying it. Every error rate previously quoted in this PR was measured under the old arrangement null, which has been deleted. Those figures — including the 0.18 familywise rate with two foci per study — described a null that no longer exists, and they do not transfer to
permute-images, whose mechanism is different. A global-null re-measurement on the current code is running; the uncorrected per-voxel rate looks near nominal in a short pilot (0.038 to 0.050), and I will post the familywise rates when the full run lands rather than state them now.Test plan
pytest nimare/tests/test_meta_effectsize.py nimare/tests/test_generate.py— 95 pass. Beyond smoke tests the suite pins the properties the model rests on: reduction to DerSimonian–Laird against PyMARE, the censored likelihood against a brute-force optimum, the null shuffling within a study and never between, a degenerate collection refused rather than given p-values, convergence alone not making a voxel significant, the threshold conversion's two designs, the prevalence held at 1 where no indicator identifies it (and the interval identity that follows), and the profile interval being unbounded wherepi = 0stands.examples/02_meta-analyses/17_plot_cbes.pyis a worked example that runs end to end in ~27s.The estimator's docstring carries the conclusions and caveats; the full measurement record — every table, every route that failed, and the symbolic proofs — lives with its scripts in a companion experiments repository rather than in the module.
🤖 Generated with Claude Code
https://claude.ai/code/session_01E37peWH4mfQi7ArYFJJSSj