Structural equation models (piecewise SEM and lavaan)

1. Binary outcome at small n (n = 40) FAILS the bias gate for the unstandardised logit path shrub -> birds (+16.9%, bias MCSE 0.03) and the unstandardised indirect fire -> birds (+15.2%). This is the known small-sample upward bias of logistic ML. Coverage was fine (0.910-0.985 across all 24 estimands), and all standardised (latent-theoretic) estimands passed (|bias| <= 9.3%). Advice: report standardised effects and their CIs, or use >= ~100 units for binary outcomes. Not fixed (Firth-type bias reduction is not implemented).
2. lavaan latent model at n = 60: one estimand (std_fertility_to_biomass) had coverage 0.900 against a 0.904 bound (bias -0.4%). Intervals for standardised lavaan parameters may be slightly narrow at small n (one borderline cell, so this is weak evidence). All raw paths, loadings and indirect effects passed.
3. Standardised count coefficients with a random intercept (observation-empirical scaling) are biased by about +7% (inside the gate). The scaling uses cor(y, mu-hat), which depends on whether mu-hat includes the random effect; the module follows piecewiseSEM (matches coefs() to 5e-5). Unstandardised paths are unbiased.
4. Unstandardised indirect effects through non-Gaussian mediators are not well-defined (products mix link scales). The module warns, and only standardised indirect effects should be compared. Count (Poisson) and binary (logit) mediators were recovery-tested at n = 150 (added 2026-10-19, both PASS 24/24): the module recovers the planted products and the standardised indirect effects, but the products themselves remain hard to interpret on the response scale.
5. Omitted paths / unmeasured confounders bias indirect effects strongly (-42% in the negative control). The d-sep test catches this omission only 30% of the time at n = 40 and 73% at n = 100.
6. Not supported: interactions or nonlinear paths, feedback loops, imperfect detection, lavaan with categorical (ordinal) endogenous variables (not validated), multigroup piecewise. Multigroup lavaan was cross-checked against separate fits (identical), but not recovery-tested.
7. The bootstrap CI recovery used 300 resamples and nsim = 100 (time). Coverage 0.90-0.99.

Status: default design PASSES the acceptance gate (all 24 estimands, 200 replicates). Two harder designs each have a documented failure (see "Known limits"). The negative control fails the gate as intended.

All runs use fast_settings (Monte Carlo CIs with 2,000 draws, d-sep on). These are likelihood fits, not MCMC, so "converged" means the optimiser reported convergence (and, for lavaan, an admissible solution); R-hat does not apply. Every replicate in every run converged (0 failures).

Reproduce: bin/ejob -c 16 -m 32 -t 4:00:00 -l logs/sem_validation2.log -- bash dev/sem/run_validation2.sh, then bin/Rscript dev/sem/dsep_conserve.R 1000 16, bash dev/sem/run_mediators.sh (count/binary mediator designs; code md5s in dev/sem/code_md5_mediators.txt) and bin/Rscript dev/sem/validation_tables.R. Outputs are in results/sem/validation/. Code version: md5s in dev/sem/code_md5_validation.txt. After these runs only output.R changed, and only in figure code (panel theme, path-diagram legend, a caption sentence). Estimates are unaffected.

Gate (core recovery_gate): |bias| <= max(10% of |truth|, 3 x bias MCSE); when truth = 0, |bias| <= max(0.1 x empirical SD, 3 x bias MCSE). Coverage of the 95% interval must be at least the binomial lower bound for nsim (0.904 at nsim = 200; 0.885 at nsim = 100), and at least 90% of replicates must converge.

Burned-forest path model, n = 100 plots: fire -> canopy (-0.6), canopy -> shrub (-0.5), moisture -> shrub (0.4), shrub -> birds (0.4, Poisson log link), fire -> birds (-0.25). Gaussian residual SD 1. Estimands: unstandardised paths (link scale), standardised paths (piecewiseSEM conventions), and direct, indirect and total effects (raw and standardised) for fire -> birds, moisture -> birds and fire -> shrub. Truth for standardised quantities comes from a 200,000-row simulation.

Summary (full table below): |relative bias| <= 2.9% for all 24 estimands; coverage 0.935-0.980. Indirect fire -> birds: truth 0.120, mean 0.1198, coverage 0.945. Indirect moisture -> birds: truth 0.160, mean 0.157, coverage 0.960. Median fit time 0.5 s.

Same design, other interval methods (both PASS 24/24):

• delta method (nsim 200): coverage 0.935-0.980; indirect fire -> birds 0.945.
• nonparametric bootstrap, 300 resamples (nsim 100): coverage 0.900-0.990 (bound 0.885).

Data generated WITH a direct fire -> shrub path (0.4) and fitted WITHOUT it (n = 100, nsim 200):

Paths not involved in the omission (fire -> canopy, shrub -> birds, fire -> birds, moisture -> shrub) still pass (|bias| <= 1.7%, coverage 0.935-0.965). The omission biases the canopy -> shrub path (canopy absorbs the fire effect) and every indirect or total effect through it.

The d-separation test detects the omission (dsep_rates.csv, 1,000 replicates each):

So Fisher's C holds about its nominal 5% size on true models (0.041-0.062) and detects an omitted path of standardised size ~0.3 in 73% of datasets at n = 100. At n = 40 it detects it in only 30%. The single claim aimed at the omitted link has more power than the global C, which dilutes it across 4 claims. "At least one claim P < 0.05" happens in 18-20% of true-model datasets, so single claims must not be over-read (the report says this).

d-sep claim direction (causal vs piecewiseSEM conserve). In all designs above piecewiseSEM's basis set contained no duplicated claims, so both options gave identical tests (rows identical in dsep_rates.csv). dev/sem/dsep_conserve.R uses a design where a claim IS duplicated (x -> m -> y1 Poisson, x -> y2 binary, y1 and y2 unlinked; 1,000 replicates; dsep_conserve.csv). True-model rejection: n = 60, causal 0.056 vs conserve 0.061; n = 150, both 0.058 (MCSE ~0.007). No meaningful size inflation from conserve here.

Default design with fire -> birds = 0 and moisture -> shrub = 0, so indirect moisture -> birds = 0 (nsim 200):

All 24 estimands pass the gate. The indirect test is conservative (0.03) because a product of two coefficients, one of which is zero, has a heavy-tailed, peaked distribution (expected for a product of coefficients).

dev/sem/crosscheck.R -> results/sem/validation/crosscheck.csv, crosscheck_ci_methods.csv, crosscheck.log.

power sem --vary n=30,50,75,100,150,200 --nsim 100 (default effects; results/sem/validation/power_n/). Power = share of 95% intervals excluding 0; coverage in brackets.

Reading: about 50 plots detect the direct and indirect fire effects, and about 100 detect the weaker moisture route. The total fire effect is small because direct and indirect partly cancel (-0.25 + 0.12), so it needs ~200 plots. Single coverages of 0.87-0.89 at nsim = 100 are within Monte Carlo error of 0.95 for one cell out of 30 (binomial 95% lower bound 0.89).

results/sem/example/sem_report.html (120 plots, Poisson GLMM birds) was read by hand. Every interpretation sentence matches the table. Examples: -0.0587 x 212 (1 SD of dNBR) = -12.4 % canopy; exp(0.0286 x 8.18) = 1.26x birds per SD of shrub; exp(0.000442 x 212) = 1.10x for the indirect route; exp(-0.000494 x 212) = 0.90x total. Fisher's C = 11.06, df 8, P = 0.198. The one flagged claim (moisture -> birds, P = 0.023) appears in text and as a dotted arrow. Conditional-effects plots are on the response scale (expected detections for the Poisson component). The lavaan example (results/sem/example_lavaan/) reports chi-square/CFI/RMSEA/SRMR, loadings, and the latent box shaded in the diagram.

bootstrap_ci (nsim = 100, failed = 0, converged = 100)

Coverage bound 0.885. Median time per fit 5.7 s.

default (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.5 s.

delta_ci (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.1 s.

glmm_poisson (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 1.0 s.

hard_small_binary (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.5 s.

lavaan_latent (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.7 s.

lavaan_latent_small (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.7 s.

misspec_omitted_path (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.5 s.

nbinom_counts (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.5 s.

null_effects (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.5 s.

Verdict: VERIFIED_WITH_LIMITS. The default design passes, and so do three new designs simulated with independent code. The controls behave as intended. The Known limits above stand.

Literature and code check. Standardisation matches piecewiseSEM 2.3.1 GetSDy/scaleGLM (binomial: latent-theoretic sqrt(var(eta) + pi^2/3), or + 1 for probit; Poisson/NB: Menard observation-empirical sd(eta)/cor(y, mu-hat); Gaussian: sd(y)). Fisher's C = -2 sum ln P on 2k df, as in fisherC (Shipley 2000, 2009). Indirect effects are sums of products over directed routes, and the Monte Carlo CIs follow Preacher & Selig (2012). simulate() and fit() share only the formula parser and the route enumerator. Because the independent runs below use none of the module's simulation or truth code, any error in that shared code would have shown up.

Builder's runs repeated (dev/sem/verify/rerun_builder.sh, results/sem/verify/builder_rerun/, nsim 200). All numbers reproduced exactly: default PASS 24/24; binary n = 40 FAIL 2/24 (+16.9%, +15.2%); omitted-path control PASS only 14/24 (indirect fire -> birds -41.6%, coverage 0.275); null effects PASS 24/24, rejection 0.055/0.045/0.030.

Independent planted truth (dev/sem/verify/verify_sem.R, which has its own simulators and truth: analytic for lavaan, a 1e6-row simulation for the non-Gaussian scalings. Output is in results/sem/verify/, log logs/sem_verify.log, nsim 200 each, 0 failed fits). The design is a grassland food chain: rain ~ N(80, 15) and grazing ~ N(2, 0.8) (non-unit scales); grass is Gaussian with a random watershed intercept fitted by lmer (15 watersheds, a component type the builder did not validate); insects is Gaussian; lizards is the sink. The effects reported are rain -> lizards (3 indirect routes, no direct path), grazing -> lizards (direct + 2 routes) and grazing -> insects.

In V1/V3/V6, path grass -> lizards (truth 0.04) averaged +15%. These runs share one set of datasets. A refit of the same datasets with MASS::glm.nb agreed to 5e-16, and 4,000 fresh glm.nb replicates gave -2.2% (MCSE 1.7%), so the +15% is Monte Carlo noise, not a module bias (dev/sem/verify/grass_bias_check.R). V5 shows again that a weak omitted path biases indirect effects strongly while the d-sep test seldom detects it.

Usability. Messy versions of the example data (dev/sem/verify/messy/) were tested.

• Missing values: check_data lists the rows.
• Wrong column name: now suggests the near match ("Did you mean shrub_cover for shrub_cover_pct?").
• Non-integer or negative counts: rows and values are now listed.
• Misspelled site IDs ('BU1', 'bu01 ' next to 'BU01'): previously NOT detected; they silently became new random-effect levels. Now flagged with row numbers.

The example report was reread and the numbers rechecked, e.g. exp(-0.000936 x 212) = 0.82, exp(0.000442 x 212) = 1.10 and exp(0.0286 x 8.18) = 1.26. The conditional-effects plots are on the response scale.

Changes made (all in output.R; estimates are unaffected):

1. check_data: near-duplicate grouping IDs (case, spaces, punctuation, leading zeros), a "did you mean" hint for missing columns, and the offending rows for bad count or binary values.

2. interpret: when direct and indirect effects have opposite signs, the text said the total "understates both". That is false when

   total	>	weaker effect	(the example:	total	0.000494 >	indirect	0.000442). It now says the total is smaller in magnitude than the stronger of the two.

3. interpret: a new "Detection versus ecology" paragraph for count or binary responses (they mix abundance or occurrence with detectability), and a note that effects in random-intercept components are for a typical group.
4. The gate description at the top of this file now matches core recovery_gate (3 MCSE, not 2; when truth = 0, max(0.1 x empirical SD, 3 MCSE)).

Remaining limits: as in Known limits 1-7. (Count and binary mediators were recovery-tested afterwards; see "Mediator designs" below.)

Run with bin/ejob -c 16 -m 32 -l logs/sem_mediators.log -- bash dev/sem/run_mediators.sh after the verification changes (code md5s in dev/sem/code_md5_mediators.txt). These close the gap noted in Known limit 4.

binary_mediator (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.4 s.

count_mediator (nsim = 200, failed = 0, converged = 200)

Coverage bound 0.904. Median time per fit 0.4 s.

Verdict: VERIFIED_WITH_LIMITS. The builder's mediator runs reproduce exactly, and two new mediator designs built with independent simulators pass. Both null controls give rejection rates of about 0.05 and the negative control fails the gate. One real interpretation problem was found and fixed (text only; estimates unchanged). Known limits 1-7 still apply.

Builder's runs repeated (dev/sem/verify/run_verify_mediators.sh, results/sem/verify/builder_rerun_mediators/, nsim 200): count_mediator PASS 24/24 and binary_mediator PASS 24/24. The summary CSVs are identical to results/sem/validation/{count,binary}_mediator/ (the only difference is an extra n_ran column written by the current core).

Independent planted truth (dev/sem/verify/verify_mediators.R; log logs/sem_verify_mediators.log; outputs results/sem/verify/mediators/; nsim 200 each; 0 failed fits, 200/200 converged). The simulators and truth are written from scratch. Standardised truth comes from a 1e6-row simulation using piecewiseSEM scalings: Gaussian sd(y); probit sqrt(var(eta) + 1); count sd(eta)/cor(y, mu). Only module$fit() is called.

The indirect tests are slightly conservative under the null (0.020-0.035), as in the builder's null_effects run.

Problem found and fixed (output.R, sem_interpret and sem_check_data).

1. With a count or binary mediator, the report text compared the raw direct, indirect and total values: "opposite signs and partly cancel, so the total (...)", and "About X% of the total runs through the mediators". The raw indirect value is a product of a coefficient on the mediator's link scale and one per observed unit of the mediator, so it is not on the same scale as the direct path. On the example data with shrub cover recoded present/absent (dev/sem/verify/check_mediator_usability.R), the raw total was +0.00073 but the standardised total was -0.20 (CI -0.37 to -0.01), and the text reported the positive raw total. For a count mediator (shrub stems), the raw indirect value was 3e-5 against a raw direct value of -9e-4, so a percentage mediated computed from raw values would be meaningless. Now, for routes through non-Gaussian mediators, the text says the raw product only tests whether the pathway is non-zero. It reports the standardised indirect and total values, and makes the cancel or percentage statements from the standardised values. For example: "The standardised direct (-0.35) and indirect (0.15) associations have opposite signs and partly cancel, so the total (-0.20)...". Text for routes through Gaussian mediators is unchanged; the example report was rechecked: exp(0.000442 x 212) = 1.10, exp(-0.000494 x 212) = 0.90, 0.0155 x 212 = 3.3.
2. If a binary mediator or response is coded as text ('yes'/'no'), check_data used to say only "response must be numeric". It now lists the values found and says how to recode them as 1/0 (in the CSV or with data: transform:).

The messy-data checks from the first pass were rerun after these edits and give the same messages: missing values (with rows), wrong column name (with a "did you mean" hint), misspelled site IDs, and non-integer or negative counts. The example fit (piecewise and lavaan) and selftest run. The conditional-effects plot for shrub -> birds is on the response scale with a 95% band; its curve ratio of 2.85 over 23-59.5% shrub matches exp(0.0286 x 36.5). Estimation code is unchanged (engine_*.R, simulate.R, graph.R, module.R have the same md5s as dev/sem/code_md5_mediators.txt), so no recovery run needed repeating. md5s after the edits are in dev/sem/verify/code_md5_after_pass2.txt.

Literature check of the mediator code. Standardised paths into a binary mediator use the latent-theoretic SD of the mediator; paths out of it use the observed SD of the 0/1 variable. This follows piecewiseSEM coefs() (Lefcheck 2016; Grace et al. 2018), and the independent truth uses the same convention. It is a convention, not a unique definition, so standardised indirect values through binary mediators are measures of relative importance, not effects on a natural scale.

Note for maintainers. While this pass was running, another agent was developing new figure code in dev/sem/figs/figures_new.R. If that code is merged into output.R, keep the sem_interpret and sem_check_data changes above (diff against dev/sem/verify/output.R.before_pass2).