bio-causal-genomics-genomic-sem

## Version Compatibility Reference examples tested with: GenomicSEM 0.0.5+ (GitHub `GenomicSEM/GenomicSEM`), lavaan 0.6-17+, LDSC v1.0.1+ (Python 3; prefer `abdenlab/ldsc-python3` v2.0.0 -- `belowlab/ldsc` v3.0.1 README states the CLI is broken; Docker `jtb114/ldsc:latest` is the belowlab fallback), baselineLD_v2.2 annotations (alkesgroup.broadinstitute.org/LDSCORE), MTAG 1.0.8+ (Python; `JonJala/mtag`), R 4.4+. Before using code patterns, verify installed versions match. If versions differ: - R: `packageVersion('GenomicSEM')` then `?ldsc`, `?commonfactor`, `?usermodel`, `?commonfactorGWAS` - Python (LDSC, MTAG): `<tool>.py -h` and inspect the source under `ldsc/` or `mtag/` GenomicSEM is GitHub-only (never on CRAN). If `ldsc()` or `usermodel()` throws an error about lavaan syntax or non-positive-definite covariance, introspect the installed API (`getMethod('ldsc')`) and adapt rather than retrying. # Genomic SEM **"Model the latent genetic architecture across several correlated GWAS"** -> Treat each GWAS as a measured indicator of one or more latent genetic factors and fit a structural equation model to the LDSC-derived genetic covariance matrix S and its sampling covariance V (Grotzinger 2019 Nat Hum Behav 3:513). The framework extends naturally to a multivariate GWAS in which a SNP is regressed on a latent factor (common-factor GWAS), with Q_SNP testing whether the SNP effect is homogeneous across factor loadings. Sample overlap between input GWAS is absorbed by the off-diagonals of V; ignoring V inflates Type-I. - R: `GenomicSEM::ldsc()` produces the (S, V) covariance pair from munged sumstats - R: `GenomicSEM::commonfactor()` fits a single-factor CFA across all traits in S - R: `GenomicSEM::usermodel()` fits an arbitrary lavaan-syntax model - R: `GenomicSEM::commonfactorGWAS()` runs SNP -> factor multivariate GWAS with Q_SNP - R: `GenomicSEM::userGWAS()` runs arbitrary multivariate SNP regression with per-path Q_SNP - Python (alternative): `mtag.py --sumstats t1,t2,t3 --out mtag_out` (multi-trait power boost on individual traits) ## Statistical Model Taxonomy | Method | Latent structure | Min traits | SNP-level test | Strength | Fails when | |--------|------------------|-----------|----------------|----------|------------| | Common-factor CFA (Grotzinger 2019) | Single F loading all traits | 3 | None (model-fit only) | Tests whether shared variance is unidimensional | Heterogeneous architecture; CFI < 0.9; near-zero loadings | | User-specified CFA (`usermodel`) | Pre-specified lavaan syntax | 3 | None | Confirmatory; arbitrary structure | Misspecified model; identification under-determined | | ESEM | Exploratory rotation; cross-loadings allowed | 6+ | None | When factor count and structure unknown | Few traits; collinear traits; rotation arbitrary | | Common-factor GWAS (`commonfactorGWAS`) | SNP -> F -> trait1..k | 3 | Wald on F + Q_SNP heterogeneity | Discovers SNPs acting via the common factor; flags Q_SNP outliers | Q_SNP-significant SNPs not interpretable as factor SNPs | | User GWAS (`userGWAS`) | Arbitrary SNP-path lavaan | 3 | Wald per path + Q_SNP | Tests SNP on any specified path | Highly parameterized models lose power | | Multivariate Wald test | Joint test across SNP -> trait paths | 2+ | Joint chi-square | Boost power when SNP affects multiple traits | Heterogeneous SNP effects collapse joint test | | Stratified GenomicSEM (Grotzinger AD et al 2022 Nat Genet 54:548) | Factor model with sLDSC-partitioned annotations | 3 | Per-annotation factor tau | Localizes heritability of the factor to functional categories | Same sLDSC failure modes (small annotation, collinearity) | | MTAG (Turley 2018 Nat Genet 50:229) | Empirical-Bayes shrinkage across correlated traits | 2 | Per-trait shrunk z-score | Boosts marginal power for any input trait | MaxFDR > 5% indicates heterogeneity violates MTAG assumption | Methodology evolves; verify the current Grotzinger 2023+ tutorials at `github.com/GenomicSEM/GenomicSEM/wiki` before locking a method. ESEM rotation choice (geomin vs target rotation) is an active area; report sensitivity to rotation. ## MTAG vs GenomicSEM Common-Factor GWAS Both methods exploit genetic correlation among input GWAS, but their goals and outputs differ. | Property | MTAG | GenomicSEM commonfactorGWAS | |----------|------|------------------------------| | Output | Per-trait shrunk z-scores | SNP effect on latent factor | | Sample-overlap handling | Bivariate LDSC intercept | Full LDSC sampling-covariance matrix V | | Heterogeneity diagnostic | MaxFDR (Turley 2018) | Q_SNP (Grotzinger 2019) | | Interpretation | "Boosted power for trait k" | "Effect on what the traits share" | | Min traits | 2 | 3 (otherwise factor not identified) | | Best when | Power-boost an individual trait | Common factor hypothesized | Both depend on accurate sampling covariance. MTAG fails (MaxFDR > 5%) under the same heterogeneity that produces large Q_SNP in GenomicSEM. The two methods should be reported together when the prior on a common factor is non-trivial; agreement increases confidence, disagreement points to architecture-specific SNPs. ## Decision Tree by Scenario | Scenario | Recommended | Why | |----------|-------------|-----| | Multi-trait GWAS power boost for one focal trait | MTAG | Optimized for per-trait marginal power | | Common-factor architecture hypothesized | GenomicSEM `commonfactorGWAS` | Tests SNP -> factor; reports Q_SNP heterogeneity | | Heterogeneous architecture (>1 latent factor) | ESEM, then confirmatory `usermodel` | Exploratory first, then confirm | | Confirming a pre-specified factor structure | `usermodel` with lavaan syntax | Confirmatory factor analysis | | Partition heritability of factor across annotations | Stratified GenomicSEM | Combines sLDSC + factor model | | Mediation in a SEM framework | `usermodel` with indirect path | Path coefficients + delta-method SE | | Sample overlap unknown or any-overlap suspected | Always use `ldsc()` output as input | V matrix off-diagonals absorb overlap | | Cross-ancestry common-factor analysis | Run per-ancestry, compare loadings; no published cross-ancestry SEM as of 2026 | Method not yet validated for mixed-ancestry V | | Single biobank for all traits (e.g., UKB only) | GenomicSEM with `ldsc()`; the V matrix will reflect overlap | Equivalent to one-sample MR -- the V matrix is the correction | | Comparing GenomicSEM and MTAG on the same traits | Run both; compare top hits + heterogeneity | Concordance increases confidence; divergence flags heterogeneity | ## Per-Method Failure Modes ### Heywood case (negative residual variance) **Trigger:** A residual variance estimate is < 0, or a standardized loading exceeds 1. **Mechanism:** Empirical underidentification; the genetic covariance matrix S is near-singular OR a trait has near-zero specific variance under the model. The maximum-likelihood / DWLS estimator runs past the boundary of the parameter space. **Symptom:** `lavaan` warning "some estimated lv variances are negative" or "covariance matrix is not positive definite"; standardized loading > 1; non-convergence. **Fix:** First, inspect the LDSC S matrix for genetic correlations near 1 (multicollinearity). Drop or merge near-identical traits. Second, constrain the offending residual variance to be non-negative in the lavaan syntax (`trait1 ~~ a*trait1; a > 0`). Third, verify the V matrix is positive definite via `chol(V_LD)`; if not, the bivariate LDSC inputs disagree on intercept signs and need re-munging. Never re-fit without diagnosing the cause. ### Sample overlap mis-specified **Trigger:** Using LDSC intercept manually or supplying covariance from non-`ldsc()` source. **Mechanism:** GenomicSEM's `ldsc()` function returns a list with `S` (genetic covariance) AND `V` (sampling covariance of the lower-triangle of S). The V off-diagonals capture sample overlap via cross-trait LDSC intercept. Skipping V and supplying only S treats all inputs as independent samples; Type-I error inflates because the sampling distribution under H0 is wrong. **Symptom:** SE on factor loadings far too small; many SNPs significant in common-factor GWAS that don't replicate; comparison to MTAG shows disagreement consistent with overlap. **Fix:** Always pass the full output of `ldsc()` -- both S and V -- to `commonfactor()`, `usermodel()`, and `commonfactorGWAS()`. Never construct S manually from rg estimates. ### Q_SNP not reported in commonfactorGWAS **Trigger:** Running `commonfactorGWAS()` and reporting only the factor p-value per SNP. **Mechanism:** Q_SNP tests heterogeneity of the SNP's effect across factor loadings (Grotzinger 2019 supplement). A SNP with significant Q_SNP violates the common-factor assumption: its effect is NOT mediated by the factor, and the factor estimate is meaningless for that SNP. **Symptom:** Top "common-factor SNPs" are dominated by trait-specific effects; replication in independent cohorts is poor for SNPs with high Q_SNP. **Fix:** Always report Q_SNP p-value alongside the factor p-value. Flag SNPs with Q_SNP p < 5e-8 / N_factor_SNPs (Bonferroni for the discovered set) as architecture-violating and exclude from "common-factor SNP" claims. Re-fit those SNPs in `userGWAS()` with separate paths to each trait. ### Trait inclusion under heterogeneous factor structure **Trigger:** Forcing a common-factor model on traits that don't share a single latent factor. **Mechanism:** When two or more traits load on a different factor than the rest, the single-factor model misfits. lavaan still returns parameter estimates but model fit is poor. **Symptom:** CFI < 0.9; RMSEA > 0.08; some standardized loadings near 0 while others near 1; chi-square highly significant even after accounting for N. **Fix:** Run ESEM first (`commonfactor` then `usermodel` with cross-loadings allowed) to discover structure. If two factors emerge, fit a two-factor `usermodel`. Drop traits with near-zero loadings on all factors. Document the model search. ### MTAG MaxFDR > 5% **Trigger:** Running MTAG on traits with low pairwise genetic correlation or with one trait that has a very different architecture. **Mechanism:** MTAG assumes a homogeneous variance-covariance structure across SNPs. When heterogeneity dominates, the empirical-Bayes shrinkage can over-claim SNPs in the focal trait. Turley 2018 defines MaxFDR as the maximum estimated false discovery rate under worst-case heterogeneity; > 5% invalidates the published trait-specific summary statistics. **Symptom:** MTAG output file reports `maxFDR` > 0.05; per-trait MTAG hits don't replicate in independent cohorts. **Fix:** Check pairwise rg via LDSC; if any pair is < 0.7, MTAG is risky. Drop the most heterogeneous trait and re-run. Alternatively, switch to GenomicSEM's `commonfactorGWAS` which models heterogeneity explicitly via Q_SNP. ### Non-positive-definite V_LD matrix **Trigger:** LDSC inputs from different ancestry GWAS, or one trait with very low mean chi-square (< 1.02). **Mechanism:** V is the sampling covariance of vech(S); when individual entries of S have huge SE relative to off-diagonal covariance, the resulting V is not positive definite (negative eigenvalues). **Symptom:** `commonfactor()` errors with "matrix is not positive definite" before fitting; `eigen(LDSCoutput$V)$values` shows negative values. **Fix:** Verify per-trait mean chi-square via `ldsc()` log; below 1.02, exclude that trait. Verify all GWAS are EUR ancestry (or match ancestry of LD scores). Apply nearest-PD smoothing via `Matrix::nearPD(V)$mat` ONLY as a last resort and document the approximation in methods. ## Model Fit Diagnostics | Index | Acceptable | Good | Source | |-------|-----------|------|--------| | CFI | >= 0.90 | >= 0.95 | Hu & Bentler 1999 Struct Equ Model 6:1 | | TLI / NNFI | >= 0.90 | >= 0.95 | Hu & Bentler 1999 | | RMSEA | <= 0.08 | <= 0.05 | Hu & Bentler 1999 | | SRMR | <= 0.08 | <= 0.05 | Hu & Bentler 1999 | | chi-square p-value | (less informative at large N) | n/a | Penalize for N inflation | AIC / BIC are used for nested-model comparison (lower is better); only compare nested models fit on the same S. ## Quantitative Thresholds | Threshold | Source | Rationale | |-----------|--------|-----------| | CFI >= 0.95 | Hu & Bentler 1999 | Conventional good fit; SEM literature default | | RMSEA <= 0.05 | Hu & Bentler 1999 | Conventional good fit | | RMSEA 0.05 - 0.08 | Hu & Bentler 1999 | Acceptable; flag as "adequate" not "good" | | Q_SNP p < 0.05 / N_SNP_factor | Grotzinger 2019 Nat Hum Behav 3:513 | Bonferroni for heterogeneity among factor-significant SNPs | | MTAG MaxFDR < 5% | Turley 2018 Nat Genet 50:229 | Above this, MTAG marginal trait results invalid | | Per-trait LDSC mean chi-square > 1.02 | LDSC documentation | Below this, V entries too noisy; factor SE inflated | | Standardized loading 0.3 - 0.9 typical | SEM conventions | < 0.3 trait loads weakly; > 0.95 may indicate over-fit / collinearity | | Min 3 traits for common factor | SEM identification | Single factor with k traits has k(k+1)/2 moments; needs k>=3 to identify | ## Standard Workflow **Goal:** Fit a common-factor model across correlated GWAS and run a multivariate GWAS on the factor with Q_SNP. **Approach:** Munge sumstats -> LDSC for (S, V) -> common-factor CFA -> inspect fit -> prepare SNPs -> common-factor GWAS -> report factor effects with Q_SNP flags. ```r library(GenomicSEM) # Step 1: Munge sumstats (one-time; produces .sumstats.gz files) files <- c('raw/trait1.txt', 'raw/trait2.txt', 'raw/trait3.txt') hm3 <- 'w_hm3.snplist' # HapMap3 SNP list trait_names <- c('trait1', 'trait2', 'trait3') N <- c(150000, 200000, 175000) munge(files = files, hm3 = hm3, trait.names = trait_names, N = N) # Step 2: LDSC produces both S (genetic covariance) and V (sampling covariance) traits <- c('trait1.sumstats.gz', 'trait2.sumstats.gz', 'trait3.sumstats.gz') ldsc_results <- ldsc( traits = traits, sample.prev = c(0.5, 0.5, NA), # case prevalence; NA for continuous population.prev = c(0.05, 0.05, NA), ld = 'eur_w_ld_chr/', wld = 'eur_w_ld_chr/', trait.names = trait_names ) # ldsc_results$S = genetic covariance; ldsc_results$V = sampling covariance # Step 3a: Common-factor CFA via DWLS cf_fit <- commonfactor(covstruc = ldsc_results, estimation = 'DWLS') print(cf_fit$modelfit) # CFI, RMSEA, SRMR, chi-square print(cf_fit$results) # loadings + SEs # Step 3b: Alternative -- user-specified two-factor model (>=3 indicators per factor) # Identification rule: each factor needs >= 3 indicators OR one anchor loading fixed # to 1 plus factor variance free. A factor with a single indicator is NOT identified. model_syntax <- ' F1 =~ NA*trait1 + trait2 + trait3 F2 =~ NA*trait4 + trait5 + trait6 F1 ~~ 1*F1 F2 ~~ 1*F2 F1 ~~ F2 ' user_fit <- usermodel(covstruc = ldsc_results, model = model_syntax, estimation = 'DWLS') ``` `estimation = 'DWLS'` (diagonally weighted least squares) is the default and is required when V is large; `'ML'` is faster but assumes a known V and can produce wrong SE under sample overlap. ### ESEM (Exploratory Factor Structure) When the factor structure is unknown, fit `usermodel()` with all loadings free across all factors, then apply a rotation post-fit. Rotation choices: **geomin oblique** (default; allows factor correlation), **target rotation** (Browne 2001 Multivariate Behav Res 36:111; uses a hypothesized loading template), **quartimin** (orthogonal; assumes factors are uncorrelated). ```r model_esem <- ' F1 =~ NA*trait1 + trait2 + trait3 + trait4 F2 =~ NA*trait1 + trait2 + trait3 + trait4 F1 ~~ 1*F1 F2 ~~ 1*F2 F1 ~~ F2 ' esem_fit <- usermodel(covstruc = ldsc_results, model = model_esem, estimation = 'DWLS') # Rotate post-fit via GPArotation::GPForth/GPFoblq or lavaan::rotate() ``` **Decision:** ESEM for K-factor exploration when structure is unknown; CFA via `usermodel()` once a structure is confirmed. Report rotation sensitivity (geomin vs target vs quartimin) and treat as exploratory. **Cross-loadings.** Cross-loadings (one trait loads on > 1 factor) are common in psychiatric and behavioral GWAS. Brown 2015 *Confirmatory Factor Analysis for Applied Research* recommends allowing cross-loadings first and using modification indices to guide simplification. Allow a cross-loading when constraining residual variance otherwise forces a Heywood case. Constrain when CFI < 0.9 and modification indices instead suggest a correlated residual between two indicators (which is the more parsimonious fix). ### userGWAS for Custom Path Models `userGWAS()` fits arbitrary lavaan-syntax SNP regressions and is the right tool when the SNP needs to be tested on multiple paths simultaneously (e.g. factor-mediated effect AND a direct effect on one indicator). ```r