Genomic SEM in umx (umxGSEM)

Psych_LDSC common-factor check against GenomicSEM

Timothy C. Bates

2026-10-11

Purpose

This vignette is a verification stub for genomic structural equation modeling (GSEM) in umx via umxGSEM().

It documents that a one-factor model on the built-in LDSC object Psych_LDSC recovers the same DWLS point estimates and standard errors (to printed precision) as GenomicSEM’s commonfactor() on the same input.

Rails, not vibes. umx builds OpenMx summary WLS with modern
mxData(numObs = N, observedStats = list(cov = S, useWeight = W, asymCov = V))
after reordering \(V\) into OpenMx residual order (var_*, poly_*_*).
Legacy OpenMx type = "acov" is not used or accepted.

Reference: Grotzinger et al. (2019). Genomic structural equation modeling. Nature Human Behaviour, 3, 513–525. doi:10.1038/s41562-019-0566-x.

Data: Psych_LDSC

Five traits from multivariate LDSC: SCZ, BIP, MDD, EA, INSOM.

data(Psych_LDSC, package = "umx")
str(Psych_LDSC, max.level = 1)
# names: V, S, I, N, m
# S is 5×5 genetic covariance; V is 15×15 sampling covariance of vech(S)
colnames(Psych_LDSC$S)
Slot Role
S Genetic covariance (heritabilities on diagonal)
V Sampling covariance of non-redundant elements of S
I LDSC intercepts
N Effective sample sizes for each S element
m SNP count used in LDSC

Fit a common factor with umxGSEM

Unit-variance identification (std.lv = TRUE by default). Path-style ~= is accepted as lavaan =~.

m1 = umxGSEM(
  model      = "g ~= SCZ + BIP + MDD + EA + INSOM",
  covstruc   = Psych_LDSC,
  estimation = "DWLS"
)
umxSummary(m1)

Optional explicit lavaan syntax (same model):

m1b = umxGSEM(
  model = "
    g =~ NA*SCZ + BIP + MDD + EA + INSOM
    g ~~ 1*g
  ",
  covstruc   = Psych_LDSC,
  estimation = "DWLS",
  std.lv     = FALSE
)

Inspect summary WLS data (modern names only — not acov):

tmx_show(m1, matrices = c("data.S", "data.V", "data.W"), report = "markdown")

Decimal match: umxGSEM vs GenomicSEM commonfactor

Verified against GenomicSEM commonfactor(Psych_LDSC, estimation = "DWLS") on the same Psych_LDSC object (unstandardized loadings and residual variances; unit variance on the factor).

Factor loadings

Path GenomicSEM umxGSEM SE (both)
g → SCZ 0.4314 0.4314 0.0408
g → BIP 0.4235 0.4235 0.0397
g → MDD 0.1310 0.1310 0.0148
g → EA 0.0087 0.0087 0.0058
g → INSOM 0.0165 0.0165 0.0055

Residual variances

Residual GenomicSEM umxGSEM
SCZ 0.0583 0.0583
BIP 0.0834 0.0834
MDD 0.0836 0.0836
EA 0.1123 0.1123
INSOM 0.0457 0.0457

Standardized loadings (GenomicSEM / same model structure) are approximately SCZ 0.87, BIP 0.83, MDD 0.41, EA 0.03, INSOM 0.08 — SCZ and BIP dominate the genetic factor; EA and INSOM barely load (as expected from S).

Fit indices (do not require decimal match)

Omnibus \(\chi^2\), CFI, and related indices can differ across packages because of how \(r'Wr\) is scaled, how the independence model is built, and dummy numObs. Parameters and SEs are the cross-check. For GSEM model comparison in umx, prefer change in CFI / SRMR and nested tests via umxCompare() (see package docs).

Package Example omnibus (illustrative)
umxGSEM \(\chi^2(5) \approx 388.5\); fit \(r'Wr \approx 407.3\) with numObs = 1
GenomicSEM commonfactor \(\chi^2(5) \approx 504\); CFI \(\approx 0.45\); SRMR \(\approx 0.16\)

Implementation notes (for maintainers)

  1. DWLS weights: \(W = \mathrm{diag}(1/\mathrm{diag}(V))\) after reordering \(V\) into OpenMx residual order (all variances, then free covariances).
  2. OpenMx data: observedStats = list(cov = S, useWeight = W, asymCov = V) — never legacy type = "acov".
  3. Identification: genetic scale prefers unit factor variance (free first loading), not unit loading.
  4. Numerics: elements of LDSC \(V\) are often \(10^{-4}\)–\(10^{-6}\); weights are large; tiny print digits can hide real structure (tmx_show auto-bumps digits for small matrices).
  5. Local minima: DWLS can sit in a secondary basin; umxGSEM defaults to tryHard = "yes" and eigen-based starts.

Re-run verification (maintainers): fit with eval = TRUE above; optionally compare to GenomicSEM commonfactor when that stack is available. No umx dependency on GenomicSEM or lavaan GSEM helpers is required for users.

SNP summary data for multivariate GWAS (pipeline roadmap)

Adding SNPs (GenomicSEM’s sumstats → commonfactorGWAS / userGWAS) needs per-trait GWAS files plus a 1000 Genomes MAF reference. Built-in Psych_LDSC / Anthro_LDSC are LDSC-only (S, V, …); they do not ship SNP betas.

How this maps to our LDSC objects

umx object Traits Closest workshop object SNP tutorial files
Psych_LDSC SCZ, BIP, MDD, EA, INSOM PRAC_COV (5-trait user model) Workshop has no bundled EA/INSOM SNP subsets for the 5-trait path model
(3-trait subset of Psych) SCZ, BIP, MDD only PSYCH_COV Yes — tutorial multivariate GWAS uses these three
Anthro_LDSC BMI, WHR, CO, Waist, Hip, Height, IHC, BL, BW Anthro_LDSC.RData in workshop materials Full anthropometric GWAS from original papers / consortia (not in the 1k SNP toy pack)

Workshop practical (IBG / Grotzinger materials):

Toy (subset) SNPs — for smoke tests only

[umxGSEM_dl_tutorial_files()] copies the IBG practical 1k-SNP files (and EUR LD scores / HapMap3) into a cache directory. Order must match ldsc / sumstats:

SCZ_subset.txt
BIP_subset.txt
MDD_subset.txt
reference.1000G.subset.txt    # drastically reduced 1000G ref for MAF

Notes from the tutorial authors:

tut = umxGSEM_dl_tutorial_files()
snps = umxGSEM_sumstats(
  files = tut$sumstats,
  ref = tut$ref,
  trait.names = c("SCZ", "BIP", "MDD"),
  se.logit = TRUE
)
gwas = umxGSEM_GWAS(covstruc = Psych_LDSC, SNPs = snps,
  traits = c("SCZ", "BIP", "MDD"), maxSNPs = 50)

Full-fat SNPs — production

  1. Download full GWAS summary statistics for each trait from the original consortia (PGC for SCZ/BIP/MDD, SSGAC for EA, etc.). Do not use “top hits only” dumps.
  2. Munge to LDSC/HapMap3 format (SNP, A1, A2, effect, P, …) with [umxGSEM_munge()]. HapMap3 list: Zenodo https://zenodo.org/records/7773502/files/w_hm3.snplist.gz (gunzip to w_hm3.snplist). Do not use [umxGSEM_dl_tutorial_files()] for this.
  3. LD scores for European-ancestry LDSC if you run ldsc locally: https://alkesgroup.broadinstitute.org/LDSCORE/eur_w_ld_chr.tar.bz2 (~500 MB → ~1.2 GB) and weights_hm3_noMHC.tgz (~150 MB). Ancestry must match. Skip if you already have a covstruc.
  4. Run multivariable LDSC → covstruc (S, V, I, N, m) — same shape as Psych_LDSC.
  5. Prepare SNP file for SEM-GWAS (sumstats-equivalent): align alleles to a 1000 Genomes Phase 3 reference used for SNP variance / MAF.
    Full ref download (GenomicSEM wiki):
    https://utexas.app.box.com/s/vkd36n197m8klbaio3yzoxsee6sxo11v
    (reference.1000G.maf.0.005.txt style; not the workshop .subset.)
  6. Trait order in munge / ldsc / sumstats / SEM-GWAS must stay identical.
  7. Binary traits: know whether SEs are logistic (se.logit = TRUE is common even when effects are ORs). Continuous traits: OLS = TRUE or standardized betas as appropriate.

Minimum columns before munge: rsID, A1 (effect allele), A2, effect (beta or OR), P. Sample size \(N\) (or effective \(N\)) is required for scaling.

umx API for SNPs (implemented)

See the toy chunk above (umxGSEM_dl_tutorial_files → umxGSEM_sumstats → umxGSEM_GWAS). Pass traits= on the LDSC object; do not V[1:6, 1:6] (wrong pairs).

Session

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: aarch64-apple-darwin23
#> Running under: macOS Golden Gate 27.2
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.6/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: Pacific/Auckland
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] umx_4.80.0           OpenMx_2.50.00.00-88
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.61           bslib_0.12.0       
#>  [4] ggplot2_4.0.3       htmlwidgets_1.6.4   visNetwork_2.1.4   
#>  [7] lattice_0.23-1      numDeriv_2016.8-1.1 quadprog_1.5-8     
#> [10] vctrs_0.7.3         tools_4.6.1         generics_0.1.4     
#> [13] stats4_4.6.1        curl_8.0.0          parallel_4.6.1     
#> [16] tibble_3.3.1        xts_0.14.3          pkgconfig_2.0.3    
#> [19] Matrix_1.7-6        RColorBrewer_1.1-3  S7_0.2.2           
#> [22] RcppParallel_6.2.1  lifecycle_1.0.5     compiler_4.6.1     
#> [25] farver_2.1.2        stringr_1.6.0       textshaping_1.0.5  
#> [28] mnormt_2.1.2        htmltools_0.5.9     sass_0.4.10        
#> [31] yaml_2.3.12         pillar_1.11.1       jquerylib_0.1.4    
#> [34] MASS_7.3-66         admisc_0.41         cachem_1.1.0       
#> [37] nlme_3.1-171        lavaan_0.7-3        polycor_0.8-2      
#> [40] tidyselect_1.2.1    digest_0.6.39       mvtnorm_1.4-2      
#> [43] stringi_1.8.9       dplyr_1.2.1         cowplot_1.2.0      
#> [46] fastmap_1.2.0       grid_4.6.1          cli_3.6.6          
#> [49] magrittr_2.0.5      DiagrammeR_1.0.12   pbivnorm_0.6.0     
#> [52] scales_1.4.0        TTR_0.24.4          rmarkdown_2.32     
#> [55] quantmod_0.4.29     otel_0.2.0          zoo_1.9-1          
#> [58] kableExtra_1.4.1    evaluate_1.0.5      knitr_1.52         
#> [61] viridisLite_0.4.3   rlang_1.3.0         Rcpp_1.1.2         
#> [64] xtable_1.8-8        glue_1.8.1          xml2_1.6.0         
#> [67] svglite_2.2.2       rstudioapi_0.19.0   jsonlite_2.0.0     
#> [70] R6_2.6.1            systemfonts_1.3.2